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

Annotation of /FigKernelScripts/connections_between.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3