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

Annotation of /FigKernelScripts/connections_between.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3