[Bio] / FigKernelScripts / svr_atomic_reg_coexp.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/svr_atomic_reg_coexp.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use Carp;
4 :     use SeedEnv;
5 :    
6 :     #
7 :     # This is a SAS Component
8 :     #
9 :    
10 :    
11 :     =head1 svr_atomic_reg_coexp
12 :    
13 :     Get functions of protein-encoding genes
14 :    
15 :     ------
16 :    
17 :     Example:
18 :    
19 :     svr_atomic_reg_coexp -g 83333.1 -ar 514
20 :    
21 :     would produce a table showing the pegs that tended to be coexpressed
22 :     with those in atomic regulon 514 of E,coli (genome 83333.1).
23 :     Each row lists a peg (which has correlation, but is not in the atomic regulon)
24 :    
25 :     The first column contains the "coexpresssed peg", the second the average
26 :     Peasrson correlation coefficient against the pegs in the atomic regulon,
27 :     and the subsequent collumns show scores against each peg in the atomic regulon.
28 :    
29 :     ------
30 :    
31 :     =head2 Command-Line Options
32 :    
33 :     =over 4
34 :    
35 :     =item -g Genome (must have expression data)
36 :    
37 :     =item -ar Atomic Regulon ID
38 :     =cut
39 :    
40 :     use SeedUtils;
41 :     use SAPserver;
42 :     my $sapO= SAPserver->new();
43 :     use Getopt::Long;
44 :     use ScriptThing;
45 :    
46 :     my $usage = "usage: svr_atomic_reg_coexp -g genome -ar AtomicRegulonID";
47 :    
48 :     my $genome;
49 :     my $ar;
50 :     my $rc = GetOptions('g=s' => \$genome,
51 :     'ar=s' => \$ar,);
52 :    
53 :     if ((! $rc) || (! $genome) || (! $ar)) { print STDERR $usage; exit }
54 :    
55 :     my $arH = $sapO->atomic_regulons( -id => $genome );
56 :     if ($arH)
57 :     {
58 :     my $fids = $arH->{"$genome:$ar"};
59 :     if ($fids)
60 :     {
61 :     my $coexpH = $sapO->coregulated_fids( -ids => $fids );
62 :     if ($coexpH)
63 :     {
64 :     my @pegs_in_ar = sort { &SeedUtils::by_fig_id($a,$b) } keys(%$coexpH);
65 :     print "\t\t". join("\t",map { $_ =~ /^\S+\.(\d+)$/; $1 } @pegs_in_ar),"\n";
66 :     my %hits;
67 :     foreach my $peg (@pegs_in_ar)
68 :     {
69 :     my $h = $coexpH->{$peg};
70 :     foreach my $peg2 (keys(%$h))
71 :     {
72 :     if ((! $coexpH->{$peg2}) && ($h->{$peg2} >= 0.5))
73 :     {
74 :     $hits{$peg2}->{$peg} = $h->{$peg2};
75 :     }
76 :     }
77 :     }
78 :     my %to_col;
79 :     for (my $i=0; ($i < @pegs_in_ar); $i++)
80 :     {
81 :     $to_col{$pegs_in_ar[$i]} = $i+2;
82 :     }
83 :    
84 :     my @rows;
85 :     foreach my $peg2 (keys(%hits))
86 :     {
87 :     my $h = $hits{$peg2};
88 :     my @in = keys(%$h);
89 :     my $row = [$peg2,0.0];
90 :     for (my $i=2; ($i < (@pegs_in_ar+2)); $i++)
91 :     {
92 :     $row->[$i] = ' ';
93 :     }
94 :     my $sc = 0;
95 :     foreach my $peg (@in)
96 :     {
97 :     $sc += $h->{$peg};
98 :     $row->[$to_col{$peg}] = sprintf("%0.3f",$h->{$peg});
99 :     }
100 :     $row->[1] = sprintf("%0.3f",$sc/@pegs_in_ar);
101 :     if ($row->[1] >= 0.5)
102 :     {
103 :     push(@rows,$row);
104 :     }
105 :     }
106 :     foreach my $row (sort { $b->[1] <=> $a->[1] } @rows)
107 :     {
108 :     print join("\t",@$row),"\n";
109 :     }
110 :     }
111 :     }
112 :     }
113 :    
114 :    
115 :    
116 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3