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

Annotation of /FigKernelScripts/check_for_new_genomes.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 :     my $chunkF = "tmp_chunk_$$.fasta";
5 :    
6 :     use FIG;
7 :     my $fig = new FIG;
8 :    
9 :     my($old,$new);
10 :     my $usage = "usage: check_new_genomes OldF NewF";
11 :     (
12 :     ($old = shift @ARGV) &&
13 :     ($new = shift @ARGV)
14 :     )
15 :     || die $usage;
16 :    
17 :     my $file_of = {};
18 :     foreach my $pair (map { $_ =~ /^(\S+)\s+(\S+)/; [$1,$2] } `cat $old $new`)
19 :     {
20 :     my($genome,$file) = @$pair;
21 :    
22 :     if ($file_of->{$genome}) { die "$genome is duplicated" }
23 :     if ($genome !~ /^\d+\.\d+$/) { die "$genome is not in the proper format" }
24 :     $file_of->{$genome} = $file;
25 :     if (! -s "$file.nsq")
26 :     {
27 :     &FIG::run("formatdb -i $file -p F");
28 :     }
29 :     }
30 :    
31 :     my $replace = {};
32 :     my @poss = &get_poss($chunkF,$new,$old);
33 :     # my @poss = map { chomp; [split(/\s+/,$_)] } `cat tmp_poss`;
34 :    
35 :     foreach my $x (@poss)
36 :     {
37 :     my($g1,$g2) = @$x;
38 :     # my $gs1 = $fig->genus_species($g1);
39 :     # my $gs2 = $fig->genus_species($g2);
40 :     # print STDERR "Checking $gs1 [$g1] and $gs2 [$g2]\n";
41 :     &test(@$x,$file_of,$replace,$fig);
42 :     }
43 :    
44 :     my %oldH = map { $_ =~ /^(\S+)/; $1 } `cat $old`;
45 :     my %newH = map { $_ =~ /^(\S+)/; $1 } `cat $new`;
46 :    
47 :     my($from,$to,$n);
48 :     foreach $_ (keys(%$replace))
49 :     {
50 :     $from = $_;
51 :     $n = 10;
52 :     while ($n-- && $replace->{$replace->{$from}})
53 :     {
54 :     $from = $replace->{$from};
55 :     }
56 :    
57 :     if ($n == 0)
58 :     {
59 :     print STDERR &Dumper(['CYCLE',$from,$replace]);
60 :     die "aborted";
61 :     }
62 :     else
63 :     {
64 :     if ($oldH{$_} && $newH{$replace->{$_}})
65 :     {
66 :     print $replace->{$_}," REPLACES $_\n";
67 :     }
68 :     else
69 :     {
70 :     print "DELETE\t$_\n";
71 :     }
72 :     }
73 :     }
74 :    
75 :     sub get_poss {
76 :     my($chunkF,$new,$old) = @_;
77 :    
78 :     my @against = map { $_ =~ /^(\S+)\s+(\S+)/; [$1,$2] } `cat $old $new`;
79 :    
80 :     my $hits = {};
81 :     foreach my $pair (map { $_ =~ /^(\S+)\s+(\S+)/; [$1,$2] } `cat $new`)
82 :     {
83 :     my($genome,$file) = @$pair;
84 :     $/ = "\n>";
85 :     open(IN,"<$file") || die "could not open $file";
86 :     my $done = 0;
87 :     while ((! $done) && (defined($_ = <IN>)))
88 :     {
89 :     chomp;
90 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
91 :     {
92 :     my $id = $1;
93 :     my $seq = $2;
94 :     $seq =~ s/\s//g;
95 :     if (length($seq) > 700)
96 :     {
97 :     $done = 1;
98 :     $seq = substr($seq,100,500);
99 :     open(CHUNK,">$chunkF") || die "could not open $chunkF";
100 :     print CHUNK ">$genome:$id\n$seq\n";
101 :     close(CHUNK);
102 :     $/ = "\n";
103 :     &process_chunk($chunkF,\@against,$hits,$genome);
104 :     unlink($chunkF);
105 :     }
106 :     $/ = "\n>";
107 :     }
108 :     }
109 :     close(IN);
110 :     $/ = "\n";
111 :     }
112 :     return sort { $a->[0] <=> $b->[0] } map { [split(/\t/,$_)] } keys(%$hits);
113 :     }
114 :    
115 :     sub process_chunk {
116 :     my($chunkF,$against,$hits,$genome) = @_;
117 :    
118 :     foreach my $pair (@$against)
119 :     {
120 :     my($genome2,$contigsF) = @$pair;
121 :     my @blastout = `blastall -m8 -i $chunkF -d $contigsF -p blastn -FF -e 1.0e-100`;
122 :     foreach my $hit (map { chomp; [split(/\t/,$_)] } @blastout)
123 :     {
124 :     my($id1,$id2,$iden,undef,undef,undef,$b1,$e1,$b2,$e2) = @$hit;
125 :     if (($genome ne $genome2) &&
126 :     ($iden >= 98) &&
127 :     (abs($e1 - $b1) > 495))
128 :     {
129 :     my($g1,$g2) = sort { $a <=> $b } ($genome,$genome2);
130 :     $hits->{"$g1\t$g2"} = 1;
131 :     }
132 :     }
133 :     }
134 :     }
135 :    
136 :     use Digest::MD5;
137 :    
138 :     sub test {
139 :     my($genome1,$genome2,$file_of,$replace,$fig) = @_;
140 :    
141 :     my ($md52id1,$id2md51) = &load_md5($file_of->{$genome1});
142 :     my ($md52id2,$id2md52) = &load_md5($file_of->{$genome2});
143 :    
144 :     my $tmpF1 = "tmp_left_1_$$.fasta";
145 :     my $tmpF2 = "tmp_left_2_$$.fasta";
146 :     my $sizes1 = &build_remaining($file_of->{$genome1},$id2md51,$md52id2,$tmpF1);
147 :     my $sizes2 = &build_remaining($file_of->{$genome2},$id2md52,$md52id1,$tmpF2);
148 :     # print STDERR &Dumper($md52id1,$id2md51,$md52id2,$id2md52,$sizes1,$sizes2);
149 :    
150 :     my($gs1,$gs2);
151 :     if ((! -s $tmpF1) && (! -s $tmpF2)) ### exactly the same
152 :     {
153 :     if ((-M $file_of->{$genome1}) <= -M $file_of->{$genome2}) ## if Genome1 is more recent
154 :     {
155 :     &replace($fig,$replace,$genome1,$genome2,1.0,1.0);
156 :     }
157 :     else
158 :     {
159 :     &replace($fig,$replace,$genome2,$genome1,1.0,1.0);
160 :     }
161 :     }
162 :     elsif (-s $tmpF1 && (! -s $tmpF2)) ### Genome1 has some extra stuff
163 :     {
164 :     &replace($fig,$replace,$genome2,$genome1,0.0,1.0);
165 :     }
166 :     elsif (-s $tmpF2 && (! -s $tmpF1)) ### Genome2 has some extra stuff
167 :     {
168 :     &replace($fig,$replace,$genome1,$genome2,1.0,0.0);
169 :     }
170 :     else
171 :     {
172 :     &resolve_with_blast($fig,$genome1,$genome2,$tmpF1,$tmpF2,$replace,$sizes1,$sizes2,$file_of);
173 :     }
174 :     system "rm $tmpF1\* $tmpF2\*";
175 :     }
176 :    
177 :     sub build_remaining {
178 :     my($file,$id2md5,$md52id,$tmpF) = @_;
179 :    
180 :     my $sizes = {};
181 :     open(IN,"<$file") || die "could not open $file";
182 :     open(OUT,">$tmpF") || die "could not open $tmpF";
183 :     $/ = "\n>";
184 :    
185 :     while (defined($_ = <IN>))
186 :     {
187 :     chomp;
188 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
189 :     {
190 :     my $id = $1;
191 :     my $seq = $2;
192 :     if (! $md52id->{$id2md5->{$id}})
193 :     {
194 :     my $ln = length($seq);
195 :     if ($ln > 200)
196 :     {
197 :     $sizes->{$id} += ($ln - 200);
198 :     print OUT ">$id\n$seq\n";
199 :     }
200 :     }
201 :     }
202 :     }
203 :     close(IN);
204 :     close(OUT);
205 :     $/ = "\n";
206 :    
207 :     return $sizes;
208 :     }
209 :    
210 :     sub resolve_with_blast {
211 :     my($fig,$genome1,$genome2,$tmpF1,$tmpF2,$replace,$sizes1,$sizes2,$file_of) = @_;
212 :    
213 :     &FIG::run("formatdb -i $tmpF1 -p F");
214 :     &FIG::run("formatdb -i $tmpF2 -p F");
215 :     my @blast1 = map { chomp; [split(/\t/,$_)] }`blastall -m8 -i $tmpF1 -d $tmpF2 -p blastn -FF -e 1.0e-100`;
216 :     my @blast2 = map { chomp; [split(/\t/,$_)] }`blastall -m8 -i $tmpF2 -d $tmpF1 -p blastn -FF -e 1.0e-100`;
217 :    
218 :     my $cov1 = &covered(\@blast1,$sizes1,$genome1);
219 :     my $cov2 = &covered(\@blast2,$sizes2,$genome2);
220 :     # print STDERR "cov1=$cov1 cov2=$cov2\n";
221 :    
222 :     my $THRESH = 0.98;
223 :     if (($cov1 >= $THRESH) && (! ($cov2 >= $THRESH))) { &replace($fig,$replace,$genome1,$genome2,$cov1,$cov2) }
224 :     elsif (($cov2 >= $THRESH) && (! ($cov1 >= $THRESH))) { &replace($fig,$replace,$genome2,$genome1,$cov2,$cov1) }
225 :     elsif (($cov1 >= $THRESH) && ($cov2 >= $THRESH))
226 :     {
227 :     my($gs1,$gs2);
228 :     if ((-M $file_of->{$genome1}) <= (-M $file_of->{$genome2}))
229 :     {
230 :     &replace($fig,$replace,$genome1,$genome2,$cov1,$cov2);
231 :     }
232 :     else
233 :     {
234 :     &replace($fig,$replace,$genome2,$genome1,$cov2,$cov1) ;
235 :     }
236 :     }
237 :     }
238 :    
239 :     sub covered {
240 :     my($blast,$sizes,$genome) = @_;
241 :    
242 :     my $all = 0;
243 :     foreach $_ (keys(%$sizes))
244 :     {
245 :     # print STDERR "contig: $_ $sizes->{$_}\n";
246 :     $all += $sizes->{$_};
247 :     }
248 :     # print STDERR &Dumper($all,$sizes);
249 :    
250 :     my @cov = sort { ($a->[0] cmp $b->[0]) or ($a->[1] <=> $b->[1]) }
251 :     map { [$_->[0],&FIG::min($_->[6],$_->[7]),&FIG::max($_->[6],$_->[7])] }
252 :     grep { $_->[2] >= 98 }
253 :     @$blast;
254 :    
255 :     my $by_contig = {};
256 :    
257 :     my($i,$j,$bound);
258 :     $i = 0;
259 :    
260 :     my $tot = 0;
261 :     while ($i < @cov)
262 :     {
263 :     # print STDERR join("\t",($cov[$i]->[0],$cov[$i]->[1],$cov[$i]->[2])),"\n";
264 :     $tot += ($cov[$i]->[2] - $cov[$i]->[1]) + 1;
265 :    
266 :     if ($i == $#cov)
267 :     {
268 :     $i++;
269 :     }
270 :     else
271 :     {
272 :     $bound = $cov[$i]->[2];
273 :     # print STDERR "bound set to $bound initially\n";
274 :     $j = $i+1;
275 :     while (($j < @cov) && ($cov[$j]->[0] eq $cov[$i]->[0]) && ($cov[$j]->[1] <= $bound))
276 :     {
277 :     if ($cov[$j]->[2] > $bound)
278 :     {
279 :     $tot += ($cov[$j]->[2] - $bound) + 1;
280 :     # print STDERR join("\t",($cov[$j]->[0],$bound,$cov[$j]->[2])),"\n";
281 :     $bound = $cov[$j]->[2];
282 :     # print STDERR "bound set to $bound\n";
283 :     }
284 :     $j++;
285 :     }
286 :     $i = $j;
287 :     }
288 :     }
289 :     my $frac = sprintf("%0.3f",$tot/$all);
290 :     print STDERR "coverage: $frac\n";
291 :     return $frac;
292 :     }
293 :    
294 :     sub load_md5 {
295 :     my($file) = @_;
296 :    
297 :     my $md52id = {};
298 :     my $id2md5 = {};
299 :    
300 :     $/ = "\n>";
301 :     open(CONTIGS,"<$file") || confess "could not open $file";
302 :     while (defined($_ = <CONTIGS>))
303 :     {
304 :     chomp;
305 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
306 :     {
307 :     my $id = $1;
308 :     my $seq = uc $2;
309 :     $seq =~ s/\s//g;
310 :     my $md5 = Digest::MD5::md5_hex( $seq );
311 :     $md52id->{$md5} = $id;
312 :     $id2md5->{$id} = $md5;
313 :     }
314 :     }
315 :     $/ = "\n";
316 :     close(CONTIGS);
317 :     return ($md52id,$id2md5);
318 :     }
319 :    
320 :     sub replace {
321 :     my($fig,$replace,$genome1,$genome2) = @_;
322 :    
323 :     # my $gs1 = $fig->genus_species($genome1);
324 :     # my $gs2 = $fig->genus_species($genome2);
325 :     # print STDERR "$genome1 [$gs1] => $genome2 [$gs2]\n";
326 :     $replace->{$genome1} = {$genome2};
327 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3