[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.1 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3