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

Annotation of /FigKernelScripts/get_coupling_values.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use SeedEnv;
3 :     use Data::Dumper;
4 : olson 1.2 use FIG_Config;
5 : overbeek 1.1
6 :     my $usage = "get_coupling_values families.2c [URL]";
7 :     my($fam2c);
8 :     (
9 :     ($fam2c = shift @ARGV)
10 :     )
11 :     || die $usage;
12 :    
13 : olson 1.2 my $genome_sets = "$FIG_Config::global/genome.sets";
14 :    
15 : overbeek 1.1 my $url = shift @ARGV;
16 :     my $sapO = new SAPserver( url => $url );
17 :    
18 :     my $tmp1 = "tmp.ff.$$.1";
19 :     my $tmp2 = "tmp.ff.$$.2";
20 :     my $tmp3 = "tmp.ff.$$.3";
21 :     &SeedUtils::run("cut -f1,2 $fam2c | sort -T . -k 2 > $tmp1");
22 :    
23 :     # my $otuMaps = $sapO->representative_genomes;
24 :     # my $to_otu = $otuMaps->[0];
25 : olson 1.2
26 :     open(GS, "<", $genome_sets) or die "Cannot open $genome_sets: $!";
27 :     my $to_otu = { map { ($_ =~ /^(\d+)\t(\d+\.\d+)/) ? ($2 => $1) : () } <GS> };
28 :     close(GS);
29 :    
30 : overbeek 1.1 my %otus = map { $to_otu->{$_} => 1 } keys(%$to_otu);
31 :    
32 :     $_ = keys(%otus);
33 :     print STDERR "$_ OTUs in analysis\n";
34 :    
35 :     open(TMP,"<",$tmp1) || die "could not open $tmp1";
36 :     my $last = <TMP>;
37 :     my $n = 0;
38 :     while (defined($last) && ($last =~ /^\S+\tfig\|(\d+\.\d+)/) && ($n < 40000))
39 :     {
40 :     my $genome = $1;
41 :     my $pegs = [];
42 :     my $to_ff = {};
43 :     while ($last && ($last =~ /^(\S+)\t(fig\|(\d+\.\d+)\.peg\.\d+)/) && ($3 eq $genome))
44 :     {
45 :     push(@$pegs,$2);
46 :     $to_ff->{$2} = $1;
47 :     $last = <TMP>;
48 :     }
49 :     my $otu = $to_otu->{$genome};
50 :    
51 :     if (defined($otu))
52 :     {
53 :     print STDERR "processing $genome\n";
54 :     &process_genome($sapO,$otu,$pegs,$to_ff,$tmp2,$tmp3);
55 :     $n++;
56 :     }
57 :     }
58 :     close(TMP);
59 :    
60 :     my $co_occH = &summarize($tmp2,4);
61 :     my $co_expH = &summarize($tmp3,2);
62 :     unlink($tmp1,$tmp2,$tmp3);
63 :    
64 :     my %pairs = map { $_ => 1 } (keys(%$co_occH),keys(%$co_expH));
65 :     foreach my $pair (sort keys(%pairs))
66 :     {
67 :     my $sc1 = $co_occH->{$pair} ? $co_occH->{$pair} : 0;
68 :     my $sc2 = $co_expH->{$pair} ? $co_expH->{$pair} : 0;
69 :     print join("\t",($pair,$sc2,$sc1)),"\n";
70 :     }
71 :    
72 :     sub summarize {
73 :     my($tmp,$min) = @_;
74 :     open(TMP,"sort $tmp |") || die "could not open $tmp";
75 :    
76 :     my $hash = {};
77 :     $last = <TMP>;
78 :     while (defined($last) && ($last =~ /^(\S+\t\S+)/))
79 :     {
80 :     my $pair = $1;
81 :     my $score = 0;
82 :     while (defined($last) && ($last =~ /^(\S+\t\S+)/) && ($1 eq $pair))
83 :     {
84 :     $score++;
85 :     $last = <TMP>;
86 :     }
87 :     if ($score >= $min)
88 :     {
89 :     $hash->{$pair} = $score;
90 :     }
91 :     }
92 :     close(TMP);
93 :     return $hash;
94 :     }
95 :    
96 :     sub process_genome {
97 :     my($sapO,$otu,$pegs,$to_ff,$tmp2,$tmp3) = @_;
98 :    
99 :     &process_contiguity($sapO,$otu,$pegs,$to_ff,$tmp2);
100 :     &process_expression($sapO,$otu,$pegs,$to_ff,$tmp3);
101 :     }
102 :    
103 :     sub process_expression {
104 :     my($sapO,$otu,$pegs,$to_ff,$tmp2) = @_;
105 :     open(TMP3,"| sort -u >> $tmp3") || die "could not open $tmp3";
106 :     my $fidH = $sapO->coregulated_fids( -ids => $pegs );
107 :     foreach my $peg (@$pegs)
108 :     {
109 :     if (my $fidHp = $fidH->{$peg})
110 :     {
111 :     my $f1 = $to_ff->{$peg};
112 :     my @hits = grep { $_->[1] >= 0.6 } map { [$_, $fidHp->{$_}] } keys(%$fidHp);
113 :     foreach my $hit (@hits)
114 :     {
115 :     my $f2 = $to_ff->{$hit->[0]};
116 :     if ($f1 && $f2)
117 :     {
118 :     if ($f1 gt $f2) { ($f1,$f2) = ($f2,$f1) }
119 :     print TMP3 join("\t",($f1,$f2,$otu)),"\n";
120 :     }
121 :     }
122 :     }
123 :     }
124 :     close(TMP3);
125 :     }
126 :    
127 :    
128 :     sub process_contiguity {
129 :     my($sapO,$otu,$pegs,$to_ff,$tmp2) = @_;
130 :     open(TMP2,"| sort -u >> $tmp2") || die "could not open $tmp2";
131 :     my $locH = $sapO->fid_locations( -ids => $pegs, -boundaries => 1 );
132 :     my @pegs_with_locs = sort { ($a->[1] cmp $b->[1]) or ($a->[2] <=> $b->[2]) }
133 :     map { my $peg = $_; my $loc = $locH->{$peg};
134 :     ($loc =~ /^\d+\.\d+:(\S+)_(\d+)([-+])(\d+)/) ? [$peg,$1,$2,$3,$4] : () } @$pegs;
135 :    
136 :     my $i;
137 :     for ($i=0; ($i < (@pegs_with_locs - 1)); $i++)
138 :     {
139 :     my $j;
140 :     for ($j=$i+1; ($j < @pegs_with_locs) && (&gap_sz($pegs_with_locs[$i],$pegs_with_locs[$j]) <= 5000); $j++)
141 :     {
142 :     my $f1 = $to_ff->{$pegs_with_locs[$i]->[0]};
143 :     my $f2 = $to_ff->{$pegs_with_locs[$j]->[0]};
144 :     ($f1 && $f2) || die "something is wrong: f1=$f1 f2=$f2";
145 :     if ($f1 gt $f2) { ($f1,$f2) = ($f2,$f1) }
146 :     print TMP2 join("\t",($f1,$f2,$otu)),"\n";
147 :     }
148 :     }
149 :     close(TMP2);
150 :     }
151 :    
152 :     sub gap_sz {
153 :     my($x,$y) = @_;
154 :    
155 :     if ($x->[1] ne $y->[1]) { return 1000000 }
156 :     my $min = ($x->[3] eq "+") ? ($x->[2] + $x->[4]) : $x->[2];
157 :     my $max = ($y->[3] eq "+") ? $y->[2] : ($y->[2] - $y->[4]);
158 :     return abs($max - $min);
159 :     }
160 :    
161 :     sub process_pair {
162 :     my($sapO,$peg1,$peg2,$to_ff) = @_;
163 :    
164 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3