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

Annotation of /FigKernelScripts/map_pegs.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 : efrank 1.1 use strict;
19 :     use FIG;
20 :    
21 :     my(@orgs,$org);
22 :     my($usage,$old_contigs,$old_tbl,$new_contigs,$new_tbl,$max);
23 :     my($nxt,$id,$pre,$rest,$old,$version,$mapF);
24 :    
25 :     $usage = "usage: map_pegs OldContigs OldTbl NewContigs NewTbl Map > adjusted.new.tbl";
26 :    
27 :     (($old_contigs = shift @ARGV) &&
28 :     ($old_tbl = shift @ARGV) &&
29 :     ($new_contigs = shift @ARGV) &&
30 :     ($new_tbl = shift @ARGV) &&
31 :     ($mapF = shift @ARGV) && open(MAP,">$mapF")
32 :     )
33 :     || die $usage;
34 :    
35 : parrello 1.3 my $fig = FIG->new();
36 :    
37 : efrank 1.1 foreach $id (map { $_ =~ /^(\S+)/; $1 ? $1 : () } `cat $old_tbl`)
38 :     {
39 :     $id =~ /peg\.(\d+)$/;
40 :     if ((! $max) || ($1 > $max)) { $max = $1 };
41 :     }
42 :     $nxt = $max+1;
43 :    
44 :     # print STDERR "nxt=$nxt\n";
45 :    
46 :     my %remap = map { $_->[0] => $_->[1] } &cluster_genes($old_contigs,$old_tbl,$new_contigs,$new_tbl);
47 :     # print STDERR &Dumper(\%remap);
48 :     foreach $_ (keys(%remap))
49 :     {
50 :     print MAP "$_\t$remap{$_}\n";
51 :     }
52 :    
53 :     foreach $_ (`cat $new_tbl`)
54 :     {
55 :     if ($_ =~ /^((fig\|\d+\.\d+\.peg\.)\d+)\t(.*)$/)
56 :     {
57 :     $id = $1;
58 :     $pre = $2;
59 :     $rest = $3;
60 :     if ($old = $remap{$id})
61 :     {
62 :     print "$old\t$rest\n";
63 :     }
64 :     else
65 :     {
66 :     print "$pre$nxt\t$rest\n";
67 :     print MAP "$id\t$pre$nxt\n";
68 :     $nxt++;
69 :     }
70 :     }
71 :     }
72 :    
73 :     sub cluster_genes {
74 :     my($old_contigs,$old_tbl,$new_contigs,$new_tbl) = @_;
75 :     my($contigs,$genes,$gene,$dna,$max,$kmer,$i,$line);
76 :     my($curr,@set,$j,$pair,$g1,$g2,@possibly_correspond,%possible);
77 :    
78 :     my $tmp_kmers = "tmp.kmers.$$";
79 :     my $clusters = [];
80 :     my @genes = ();
81 :    
82 :     open(KMERS,"| sort +0 -1 +1n -2 > $tmp_kmers") || die "aborted";
83 :     foreach $version (([$new_contigs,$new_tbl,"new"],[$old_contigs,$old_tbl,"old"]))
84 :     {
85 :     my($contigsF,$tblF,$which) = @$version;
86 :     # print STDERR "loading $which\n";
87 :    
88 :     if (($contigs = &get_contigs($contigsF)) &&
89 :     ($genes = &get_genes($contigs,$tblF,$which)))
90 :     {
91 :     # $_ = @$genes; print STDERR "$_ genes loaded\n";
92 :    
93 :     foreach $gene (@$genes) # $gene is [Id,Contig,Beg,End,DNA,which]
94 :     {
95 :     $dna = $gene->[4];
96 :     push(@genes,$gene);
97 :     $max = length($dna) - 20;
98 :     for ($i=0; ($i <= $max); $i++)
99 :     {
100 :     $kmer = uc substr($dna,$i,20);
101 :     $kmer =~ tr/U/T/;
102 :     if ($kmer !~ /[^ACGT]/)
103 :     {
104 :     print KMERS "$kmer\t$#genes\n";
105 :     }
106 :     }
107 :     }
108 :     }
109 :     # print STDERR "finishished loading $which\n";
110 :     }
111 :     close(KMERS);
112 :    
113 :     # print STDERR "processing kmers\n";
114 :     open(KMERS,"<$tmp_kmers")
115 :     || die "could not open $tmp_kmers";
116 :    
117 :     $line = <KMERS>;
118 :     while ($line && ($line =~ /^(\S+)\t(\d+)$/))
119 :     {
120 :     $curr = $1;
121 :     @set = ();
122 :     while ($line && ($line =~ /^$curr\t(\d+)/))
123 :     {
124 :     push(@set,$1);
125 :     $line = <KMERS>;
126 :     }
127 :    
128 :     for ($i=0; ($i < $#set); $i++)
129 :     {
130 :     for ($j=$i+1; ($j < @set); $j++)
131 :     {
132 :     if ($genes[$set[$i]]->[5] ne $genes[$set[$j]]->[5])
133 :     {
134 :     $possible{"$set[$i],$set[$j]"}++;
135 :     }
136 :     }
137 :     }
138 :     }
139 :     close(KMERS);
140 :     unlink($tmp_kmers);
141 :     # print STDERR "built kmer matches\n";
142 :    
143 :     foreach $pair (keys(%possible))
144 :     {
145 :     if ($possible{$pair} > 4)
146 :     {
147 :     ($g1,$g2) = split(/,/,$pair);
148 :     # print STDERR "possibly related: $genes[$g1]->[0] $genes[$g2]->[0]\n";
149 :     if (&possibly_correspond($genes[$g1]->[4],$genes[$g2]->[4]))
150 :     {
151 :     push(@possibly_correspond,[$g1,$g2]);
152 :     }
153 :     }
154 :     }
155 :     return map { [$genes[$_->[0]]->[0],$genes[$_->[1]]->[0]] } &compute_correspondence(\@possibly_correspond,\@genes);
156 :     }
157 :    
158 :     sub compute_correspondence {
159 :     my($possible,$genes) = @_;
160 :     my($clearly,$possibly,$pair);
161 :    
162 :     ($clearly,$possibly) = &map_unique($possible);
163 :     # print STDERR &Dumper(["possibly",$possibly]);
164 :    
165 :     # print STDERR "CLEARLY\n";
166 :     foreach $pair (@$clearly)
167 :     {
168 :     my($g1,$g2) = @$pair;
169 :     # print STDERR &Dumper($genes->[$g1],$genes->[$g2]);
170 :     }
171 :     # print STDERR "============\n";
172 :    
173 :     foreach $pair (@$possibly)
174 :     {
175 :     # print STDERR "checking for possibility: ",join(" ",@$pair),"\n";
176 :     if (&using_bounds($pair,$possibly,$clearly,$genes))
177 :     {
178 :     # print STDERR "adding possible pair to clearly ",join(" ",@$pair),"\n";
179 :     push(@$clearly,$pair);
180 :     }
181 :     else
182 :     {
183 :     # print STDERR "failed: ",&Dumper($genes->[$pair->[0]],$genes->[$pair->[1]]),"\n";
184 :     }
185 :     }
186 :     return @$clearly;
187 :     }
188 :    
189 :     sub using_bounds {
190 :     my($pair,$possibly,$clearly,$genes) = @_;
191 :     my($g1,$g2,$g1A,$g2A,$pairA,$lg1,$rg1,$lg2,$rg2);
192 :    
193 :     # print STDERR "USING BOUNDS\n";
194 :     ($g1,$g2) = @$pair;
195 :     # print &Dumper($genes->[$g1],$genes->[$g2]);
196 :     foreach $pairA (@$clearly)
197 :     {
198 :     ($g1A,$g2A) = @$pairA;
199 :     if (&better_left($genes,$g1,$lg1,$g1A) && &better_left($genes,$g2,$lg2,$g2A))
200 :     {
201 :     $lg1 = $g1A;
202 :     $lg2 = $g2A;
203 :     }
204 :    
205 :     if (&better_right($genes,$g1,$rg1,$g1A) && &better_right($genes,$g2,$rg2,$g2A))
206 :     {
207 :     $rg1 = $g1A;
208 :     $rg2 = $g2A;
209 :     }
210 :     }
211 :     # print STDERR &Dumper($lg1,$g1,$rg1);
212 :     if (defined($lg1) && defined($g1) && defined($rg1))
213 :     {
214 :     # print STDERR &Dumper($genes->[$lg1],$genes->[$g1],$genes->[$rg1]);
215 :     }
216 :    
217 :     # print STDERR &Dumper($lg2,$g2,$rg2);
218 :     if (defined($lg2) && defined($g2) && defined($rg2))
219 :     {
220 :     # print STDERR &Dumper($genes->[$lg2],$genes->[$g2],$genes->[$rg2]);
221 :     }
222 :    
223 :     if (defined($lg1) && defined($rg1) && defined($lg2) && defined($rg2))
224 :     {
225 :     foreach $pairA (@$possibly)
226 :     {
227 :     ($g1A,$g2A) = @$pairA;
228 :     if (($g1A eq $g1) && ($g2A ne $g2) && &in_bounds($genes,$g2A,$lg2,$rg2)) { return 0 }
229 :     if (($g2A eq $g2) && ($g1A ne $g1) && &in_bounds($genes,$g1A,$lg1,$rg1)) { return 0 }
230 :     }
231 :     # print STDERR "OK: $g1 $g2\n";
232 :     return 1;
233 :     }
234 :     else
235 :     {
236 :     # print STDERR "$g1 $g2 FAILED\n";
237 :     return 0;
238 :     }
239 :     }
240 :    
241 :     sub in_bounds {
242 :     my($genes,$g,$lg,$rg) = @_;
243 :    
244 :     return (($genes->[$g]->[1] eq $genes->[$lg]->[1]) &&
245 :     &left($genes->[$lg]->[2],$genes->[$lg]->[3],$genes->[$g]->[2],$genes->[$g]->[3]) &&
246 :     &left($genes->[$g]->[2],$genes->[$g]->[3],$genes->[$rg]->[2],$genes->[$rg]->[3]));
247 :     }
248 :    
249 :     sub better_left {
250 :     my($genes,$g,$lg,$ga) = @_;
251 :     my($rc);
252 :    
253 :     if (($genes->[$g]->[1] ne $genes->[$ga]->[1]) ||
254 :     (! &left($genes->[$ga]->[2],$genes->[$ga]->[3],$genes->[$g]->[2],$genes->[$g]->[3])))
255 :     {
256 :     $rc = 0;
257 :     }
258 :     elsif ((! $lg) || &left($genes->[$lg]->[2],$genes->[$lg]->[3],$genes->[$ga]->[2],$genes->[$ga]->[3]))
259 :     {
260 :     $rc = 1;
261 :     }
262 :     else
263 :     {
264 :     $rc = 0;
265 :     }
266 :     # print STDERR "better left: $g $ga = $rc\n";
267 :     return $rc;
268 :     }
269 :    
270 :     sub left {
271 :     my($b1,$e1,$b2,$e2) = @_;
272 :     my($min,$max);
273 :    
274 :     return (&FIG::max($b1,$e1) < &FIG::min($b2,$e2));
275 :     }
276 :    
277 :     sub better_right {
278 :     my($genes,$g,$rg,$ga) = @_;
279 :    
280 :     if (($genes->[$g]->[1] ne $genes->[$ga]->[1]) ||
281 :     (! &left($genes->[$g]->[2],$genes->[$g]->[3],$genes->[$ga]->[2],$genes->[$ga]->[3]))) { return 0 }
282 :     if ((! $rg) || &left($genes->[$ga]->[2],$genes->[$ga]->[3],$genes->[$rg]->[2],$genes->[$rg]->[3])) { return 1 }
283 :     return 0;
284 :     }
285 :    
286 :     sub map_unique {
287 :     my($possible) = @_;
288 :     my($pair,%hits1,%hits2,$g1,$g2);
289 :     my $clearly = [];
290 :     my $possibly = [];
291 :    
292 :     foreach $pair (@$possible)
293 :     {
294 :     ($g1,$g2) = @$pair;
295 :     $hits1{$g1}++;
296 :     $hits2{$g2}++;
297 :     }
298 :    
299 :     foreach $pair (@$possible)
300 :     {
301 :     ($g1,$g2) = @$pair;
302 :     if (($hits1{$g1} == 1) && ($hits2{$g2} == 1))
303 :     {
304 :     push(@$clearly,$pair);
305 :     }
306 :     else
307 :     {
308 :     push(@$possibly,$pair);
309 :     }
310 :     }
311 :     return ($clearly,$possibly);
312 :     }
313 :    
314 :     sub possibly_correspond {
315 :     my($dna1,$dna2) = @_;
316 :     my(@ambig,$best,$where1,$where2,$i1,$i2,$i,$max,$chunk);
317 :    
318 :     if (($dna1 =~ /[^ACGT]{4}/) || ($dna2 =~ /[^ACGT]{4}/)) { return 0 }
319 :    
320 :     if (length($dna1) > length($dna2))
321 :     {
322 :     ($dna1,$dna2) = ($dna2,$dna1);
323 :     }
324 :    
325 :     @ambig = (0);
326 :     while ($dna1 =~ /[MRWSYKBDHVN]/g)
327 :     {
328 :     push(@ambig,pos($dna1));
329 :     }
330 :     push(@ambig,length($dna1));
331 :    
332 :     $best = 0;
333 :     for ($i=0; ($i < $#ambig); $i++)
334 :     {
335 :     if (($ambig[$i+1] - $ambig[$i]) > $best)
336 :     {
337 :     $best = $ambig[$i+1] - $ambig[$i];
338 :     $where1 = $ambig[$i];
339 :     }
340 :     }
341 :     $chunk = substr($dna1,$where1,$best);
342 :     $chunk =~ s/[MRWSYKBDHVN]/./g;
343 :    
344 :     if ($dna2 =~ /$chunk/g)
345 :     {
346 :     $where2 = pos($dna2) - length($chunk);
347 :     $max = length($dna1);
348 :     if ($where2 >= $where1)
349 :     {
350 :     for ($i1=0, $i2=$where2-$where1; ($i1 < $max) && &match(substr($dna1,$i1,1),substr($dna2,$i2,1)); $i1++, $i2++) {}
351 :     return ($i1 == $max);
352 :     }
353 :     }
354 :     return 0;
355 :     }
356 :    
357 :     sub match {
358 :     my($c1,$c2) = @_;
359 :    
360 :     if ($c1 eq $c2) { return 1 }
361 :     if ($c1 !~ /^[ACGT]/) { ($c1,$c2) = ($c2,$c1) }
362 :     if ($c1 !~ /^[ACGT]/) { return 0 }
363 :     if ($c1 eq "A") { return ($c2 =~ /[MRWDHVN]/) }
364 :     elsif ($c1 eq "C") { return ($c2 =~ /[MSYBHVN]/) }
365 :     elsif ($c1 eq "G") { return ($c2 =~ /[RSKBDVN]/) }
366 :     elsif ($c1 eq "T") { return ($c2 =~ /[WYKBDHN]/) }
367 :     return 0;
368 :     }
369 :    
370 :     sub get_contigs {
371 :     my($file) = @_;
372 :     my($contigs,$id,$seq);
373 :    
374 :     open(CONTIGS,"<$file") || die "could not open $file";
375 :     print STDERR "loading contigs from $file\n";
376 :    
377 :     $contigs = {};
378 :     $/ = "\n>";
379 :     while (defined($_ = <CONTIGS>))
380 :     {
381 :     chomp;
382 :     if ($_ =~ /^>?(\S+)[^\n]*\n(.*)/s)
383 :     {
384 :     $id = $1;
385 :     $seq = $2;
386 :     $seq =~ s/\s//gs;
387 :     $contigs->{$id} = $seq;
388 :     }
389 :     }
390 :     $/ = "\n";
391 :     close(CONTIGS);
392 :     return $contigs;
393 :     }
394 :    
395 :     sub get_genes {
396 :     my($contigs,$tblF,$which) = @_;
397 :     my($contig,$loc,$dna);
398 :     my($beg,$end,$line,$seq);
399 :    
400 :     open(TBL,"<$tblF") || die "could not open $tblF";
401 :     my $genes = [];
402 :     while (defined($line = <TBL>))
403 :     {
404 :     if ($line =~ /^(\S+)\t(\S+)/)
405 :     {
406 :     $id = $1;
407 :     $loc = $2;
408 : parrello 1.3 ($contig,$beg,$end) = $fig->boundaries_of($loc);
409 : efrank 1.1 $seq = $contigs->{$contig};
410 :     if ($beg < $end)
411 :     {
412 :     if (($beg > 0) && ($end <= length($seq)))
413 :     {
414 :     $dna = uc substr($seq,$beg-1,($end+1)-$beg);
415 :     push(@$genes,[$id,$contig,$beg,$end,$dna,$which]);
416 :     }
417 :     }
418 :     else
419 :     {
420 :     if (($end > 0) && ($beg <= length($seq)))
421 :     {
422 :     $dna = uc &reverse_complement(substr($seq,$end-1,($beg+1)-$end));
423 :     push(@$genes,[$id,$contig,$beg,$end,$dna,$which]);
424 :     }
425 :     }
426 :     }
427 :     }
428 :     return $genes;
429 :     }
430 :    
431 :     sub reverse_complement {
432 :     my($seq) = (@_ == 1) ? $_[0] : $_[1];
433 :     my($rev);
434 :    
435 :     $rev = reverse($seq);
436 :     $rev =~ tr/a-z/A-Z/;
437 :     $rev =~ tr/ACGTUMRWSYKBDHV/TGCAAKYWSRMVHDB/;
438 :     return $rev;
439 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3