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

Annotation of /FigKernelScripts/connections_between.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.8 ########################################################################
2 : olson 1.2 #
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     #
18 :    
19 : overbeek 1.1
20 :     use strict;
21 :    
22 :     use FIG;
23 : olson 1.6 use FIGV;
24 :    
25 :     my $fig;
26 : overbeek 1.1
27 : overbeek 1.8 my($simsD,$iden_diff);
28 : overbeek 1.1
29 : overbeek 1.8
30 :     my $usage = "usage: connections_between [-orgdir orgdir] [-cluster N1] [-insub N2] Sims IdenDiff > connections";
31 : olson 1.6
32 :     my $orgdir;
33 : overbeek 1.8 my $must_cluster = 0;
34 :     my $must_be_in_sub = 0;
35 :    
36 : olson 1.6 while ((@ARGV > 0) && ($ARGV[0] =~ /^-/))
37 :     {
38 :     my $arg = shift @ARGV;
39 :     if ($arg =~ /^-orgdir/)
40 :     {
41 :     $orgdir = shift @ARGV;
42 :     }
43 : overbeek 1.8 elsif ($arg =~ /-cluster/)
44 :     {
45 :     $must_cluster = shift @ARGV;
46 :     }
47 :     elsif ($arg =~ /-insub/)
48 :     {
49 :     $must_be_in_sub = shift @ARGV;
50 :     }
51 : olson 1.6 else
52 :     {
53 :     die "Unknown option $arg\n";
54 :     }
55 :     }
56 :    
57 :     if ($orgdir)
58 :     {
59 :     $fig = new FIGV($orgdir);
60 :     }
61 :     else
62 :     {
63 :     $fig = new FIG;
64 :     }
65 :    
66 : overbeek 1.1 (
67 :     ($simsD = shift @ARGV) &&
68 : overbeek 1.8 ($iden_diff = shift @ARGV)
69 : overbeek 1.1 )
70 :     || die $usage;
71 :    
72 :     #
73 :     # Here we compute whether G1 <-> G2 by first finding out
74 :     #
75 :     # 1. Are they BBHs with a safety threshold of IdenDiff (using identify, rather than Psc)?
76 :     # If they are, then they are matched.
77 :     #
78 :     # 2. If
79 :     # a) G1 is not a safe BBH of G2 (but is similar)
80 :     # b) there is no other safe BBH to either
81 :     # c) they have a superior context (by at least NeighDiff)
82 :     #
83 :     # then we call it a match
84 :     #
85 :     # We compute all binary matches.
86 :     #
87 :     # Then, we compute families as transitive close over the binary matches.
88 :     #
89 :     # Any resulting family with more than one gene from a single genome is registered
90 :     # as an error.
91 :    
92 :     opendir(SIMS,$simsD) || die "Where are the similarities?";
93 :     my @sims_files = grep { $_ !~ /^\./ } readdir(SIMS);
94 :     closedir(SIMS);
95 :    
96 : overbeek 1.7 my($file,%genomes,%sims,%paralogs);
97 :     foreach $file (@sims_files)
98 :     {
99 :     if ($file =~ /^(\d+\.\d+)-(\d+\.\d+)$/)
100 :     {
101 :     my $genome1 = $1;
102 :     my $genome2 = $2;
103 :     if ($genome1 eq $genome2)
104 :     {
105 :     open(PAR,"cut -f1 $simsD/$file |") || die "could not open $simsD/$file";
106 :     while (defined($_ = <PAR>))
107 :     {
108 :     chomp;
109 :     $paralogs{$_} = 1;
110 :     }
111 :     close(PAR);
112 :     }
113 :     }
114 :     }
115 :    
116 : overbeek 1.1 foreach $file (@sims_files)
117 :     {
118 :     if ($file =~ /^(\d+\.\d+)-(\d+\.\d+)$/)
119 :     {
120 :     my $genome1 = $1;
121 :     my $genome2 = $2;
122 :     $genomes{$genome1} = 1;
123 :     $genomes{$genome2} = 1;
124 :    
125 :     if (open(SIMS,"<$simsD/$file"))
126 :     {
127 :     while (defined($_ = <SIMS>))
128 :     {
129 : overbeek 1.3 chomp;
130 :     my($peg1,$peg2,$iden,$b1,$e1,$b2,$e2,$sc,$ln1,$ln2) = split(/\t/,$_);
131 : overbeek 1.7 if (((($e1-$b1)/$ln1) > 0.7) && ((($e2-$b2)/$ln2) > 0.8) && (! $paralogs{$peg1}) && (! $paralogs{$peg2}))
132 : overbeek 1.1 {
133 : overbeek 1.3 push(@{$sims{"$genome1,$genome2"}->{$peg1}},[$iden,$peg2]);
134 : overbeek 1.1 }
135 :     }
136 :     close(SIMS);
137 :     }
138 :     else
139 :     {
140 :     die "could not open $simsD/$file";
141 :     }
142 :     }
143 :     else
144 :     {
145 :     die "$simsD/$file is not a valid similarity file";
146 :     }
147 :     }
148 :    
149 :     my @genomes = sort { $a <=> $b } keys(%genomes);
150 :     my $n = @genomes;
151 :     print STDERR "Building families from $n genomes\n";
152 :     my $genome;
153 :     foreach $genome (@genomes)
154 :     {
155 :     my $gs = $fig->genus_species($genome);
156 :     print STDERR "$genome\t$gs\n";
157 :     }
158 :     print STDERR "\n";
159 :    
160 :     my $tmpF = "/tmp/connections$$";
161 :     open(TMP,">$tmpF") || die "could not open $tmpF";
162 :     my $neighbors = {};
163 :    
164 :     my($genome1,$genome2);
165 :     foreach $genome1 (@genomes)
166 :     {
167 :     foreach $genome2 (@genomes)
168 :     {
169 : overbeek 1.5 if ( $genome1 lt $genome2 )
170 : overbeek 1.1 {
171 :     my $genomes = "$genome1,$genome2";
172 :     print STDERR "processing $genomes\n";
173 :     my $bbhs = &bbhs($sims{$genomes},$sims{"$genome2,$genome1"},$iden_diff);
174 :    
175 :     my($peg1,$peg2,);
176 :    
177 :     foreach $peg1 (sort { &FIG::by_fig_id($a,$b) } keys(%{$sims{$genomes}}))
178 :     {
179 :     my $conn;
180 :     if ($peg2 = $bbhs->{$peg1})
181 :     {
182 : overbeek 1.8 if (! $must_cluster)
183 : overbeek 1.1 {
184 : overbeek 1.8 print TMP "$peg1\t$peg2\n";
185 :     }
186 :     else
187 :     {
188 :     my $clustered = &clustered($peg1,$peg2,\%sims,$genome1,$genome2,$bbhs);
189 :     if (@$clustered >= $must_cluster)
190 : overbeek 1.1 {
191 : overbeek 1.8 my $common_subs;
192 :     if ((! $must_be_in_sub) || (@{($common_subs = &insub($fig,$clustered))} >= $must_be_in_sub))
193 :     {
194 :     print TMP "$peg1\t$peg2\n";
195 :     print STDERR "connection\t$peg1\t$peg2\t";
196 :     foreach my $x (@$clustered) { print STDERR join(",",@$x),";" } print STDERR "\t";
197 :     foreach my $x (@$common_subs) { print STDERR join(",",@$x),";" } print STDERR "\n";
198 :     }
199 : overbeek 1.1 }
200 :     }
201 :     }
202 :     }
203 :     }
204 :     }
205 :     }
206 :     close(TMP);
207 :     open(CLUSTERS,"cluster_objects < $tmpF |") || die "could not run cluster_objects on $tmpF";
208 :     while (defined($_ = <CLUSTERS>))
209 :     {
210 :     chop;
211 :    
212 :     my($peg,%genomes);
213 :     my @pegs = split(/\t/,$_);
214 :    
215 :     undef %genomes;
216 :     my $ok = 1;
217 :     foreach $peg (@pegs)
218 :     {
219 :     if ($genomes{&FIG::genome_of($peg)}) { $ok = 0 }
220 :     $genomes{&FIG::genome_of($peg)} = 1;
221 :     }
222 :    
223 :     if ($ok)
224 :     {
225 :     print "$_\n";
226 :     }
227 :     else
228 :     {
229 :     print STDERR "$_\n";
230 :     }
231 :     }
232 : overbeek 1.8 unlink($tmpF);
233 : overbeek 1.1
234 : overbeek 1.8 sub clustered {
235 :     my($peg1,$peg2,$sims,$genome1,$genome2,$bbhs) = @_;
236 :     my $neigh1 = &neighbors_of($peg1,$neighbors);
237 :     my $neigh2 = &neighbors_of($peg2,$neighbors);
238 :     my %neigh2H = map { $_ => 1 } @$neigh2;
239 :     my $matches = [];
240 :     my ($peg_in1,$peg_in2);
241 :     foreach $peg_in1 (@$neigh1)
242 : overbeek 1.1 {
243 : overbeek 1.8 $peg_in2 = $bbhs->{$peg_in1};
244 :     if ($neigh2H{$peg_in2})
245 : overbeek 1.1 {
246 : overbeek 1.8 push(@$matches,[$peg_in1,$peg_in2]);
247 : overbeek 1.1 }
248 :     }
249 : overbeek 1.8 return $matches;
250 :     }
251 : overbeek 1.1
252 : overbeek 1.8 sub insub {
253 :     my($fig,$clustered) = @_;
254 :    
255 :     my $n = [];
256 :     my $tuple;
257 :     foreach $tuple (@$clustered)
258 : overbeek 1.1 {
259 : overbeek 1.8 my($peg1,$peg2) = @$tuple;
260 :     my $func1 = $fig->function_of($peg1);
261 :     my $func2 = $fig->function_of($peg2);
262 :     if ($func1 eq $func2)
263 : overbeek 1.1 {
264 : overbeek 1.8 my @tmp = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg1);
265 :     if (@tmp > 0)
266 :     {
267 :     @tmp = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg2);
268 :     if (@tmp > 0)
269 :     {
270 :     push(@$n,[$peg1,$peg2,$func1]);
271 :     }
272 :     }
273 : overbeek 1.1 }
274 :     }
275 :     return $n;
276 :     }
277 : overbeek 1.8
278 : overbeek 1.1
279 :     sub neighbors_of {
280 :     my($peg,$neighbors) = @_;
281 :    
282 :     if (! ($_ = $neighbors->{$peg}))
283 :     {
284 :     $neighbors->{$peg} = [grep { $_ =~ /peg/ } $fig->close_genes($peg,5000)];
285 :     }
286 :     return $neighbors->{$peg};
287 :     }
288 :    
289 :     sub bbhs {
290 :     my($sims1,$sims2,$iden_diff) = @_;
291 :     my($peg1,$best1,$best2);
292 :    
293 :     my $bbhs = {};
294 :     foreach $peg1 (sort { &FIG::by_fig_id($a,$b) } keys(%$sims1))
295 :     {
296 :     if (($best1 = &best($sims1->{$peg1},$iden_diff)) &&
297 :     ($best2 = &best($sims2->{$best1},$iden_diff)) &&
298 :     ($best2 eq $peg1))
299 :     {
300 :     $bbhs->{$peg1} = $best1;
301 :     $bbhs->{$best1} = $peg1;
302 :     }
303 :     }
304 :     return $bbhs;
305 :     }
306 :    
307 :     sub best {
308 :     my($sims,$iden_diff) = @_;
309 :    
310 :     if (! $sims) { return undef }
311 :    
312 :     if ((@$sims == 1) ||
313 :     ((@$sims > 1) && (($sims->[0]->[0] - $sims->[1]->[0]) >= $iden_diff)))
314 :     {
315 :     return $sims->[0]->[1];
316 :     }
317 :     return undef;
318 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3