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

Diff of /FigKernelScripts/build_initial_objects_for_start_predictions.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.3, Thu May 12 20:38:37 2005 UTC revision 1.12, Sun Jan 20 21:44:52 2008 UTC
# Line 1  Line 1 
1    # -*- perl -*-
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  use FIG;  use FIG;
20  my $fig = new FIG;  my $fig = new FIG;
21    use FIGV;
22    
23  $usage = "usage: build_initial_objects_for_start_predictions < pegs_in_cluster > initial_objects";  $0 =~ m/([^\/]+)$/;  $self = $1;
24    warn "$0 " . join(" ", @ARGV) . "\n";
25    $usage = "build_initial_objects_for_start_predictions [-orgdir=OrgDir] [-genetic_code=CodeNum] [tbl] < pegs_to_call > initial_objects";
26    
27    use constant FID    =>  0;
28    use constant LOCUS  =>  1;
29    use constant CONTIG =>  2;
30    use constant START  =>  3;
31    use constant STOP   =>  4;
32    use constant LEFT   =>  5;
33    use constant RIGHT  =>  6;
34    use constant LEN    =>  7;
35    use constant STRAND =>  8;
36    
37    $tbl_file = "";
38    $trouble  = 0;
39    $genetic_code_number = 11;
40    while (@ARGV)
41    {
42        if ($ARGV[0] =~ m/-h(elp)?/) {
43            die "\n\tusage: $usage\n\n";
44        }
45        elsif ($ARGV[0] =~ /-orgdir=(\S+)/) {
46            warn "Using $ARGV[0]\n";
47            $fig = new FIGV($1);
48        }
49        elsif ($ARGV[0] =~ /-genetic_code=(\d+)/) {
50            warn "$self is using $ARGV[0]\n";
51            $genetic_code_number = $1;
52        }
53        elsif (-s $ARGV[0]) {
54            $tbl_file = $ARGV[0];
55        }
56        else {
57            warn "Invalid argument $ARGV[0]\n";
58        }
59        shift @ARGV;
60    }
61    die "\nusage: $usage\n\n" if $trouble;
62    
63    $code_map = &FIG::genetic_code($genetic_code_number);
64    
65  @pegs = <STDIN>;  @pegs = <STDIN>;
66  chomp @pegs;  chomp @pegs;
67  # print "NPEGS=", scalar(@pegs), "\n";  # print "NPEGS=", scalar(@pegs), "\n";
68  for ($i=0; $i < @pegs; $i++)  # strip off all but peg ids at front  
69    if ($tbl_file)
70    {
71        open(TBL, "<$tbl_file") || die "Could not read-open $tbl_file";
72        while (defined($entry = <TBL>))
73  {  {
74      if ($pegs[$i] =~ /(\S+)\t.*/)          if ($entry =~ m/^(\S+)\t(\S+)/)
75      {      {
76          $pegs[$i] = $1;              ($fid, $loc) = ($1, $2);
77                ($contig, $beg, $end, $left, $right, $len, $strand) = &from_locus($loc);
78                $tbl->{$contig}->{$strand.$end}
79                    = [ $fid, $loc, $contig, $beg, $end, $left, $right, $len, $strand ];
80            }
81            else
82            {
83                print STDERR "Skipping invalid tbl entry: $entry";
84            }
85      }      }
86  }  }
87    
88    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
89    #...If there are "new" PEGs, overwrite any existing tbl entries...
90    #-----------------------------------------------------------------------
91  foreach $peg (@pegs)  foreach $peg (@pegs)
92  {  {
93        if ($peg =~ /^(new\|(\d+\.\d+)\.peg\.\d+)\t(\S+)/)
94        {
95            $fid    = $1;
96            $genome = $2;
97            $loc    = $3;
98    
99            ($contig, $beg, $end, $left, $right, $len, $strand) = &from_locus($loc);
100            $tbl->{$contig}->{$strand.$end}
101                = [ $fid, $loc, $contig, $beg, $end, $left, $right, $len, $strand ];
102        }
103    }
104    
105    foreach $contig (keys %$tbl)
106    {
107        @$x = sort { $a->[LEFT] <=> $b->[LEFT] } values % { $tbl->{$contig} };
108    
109        for ($i=0; $i < @$x; ++$i)
110        {
111            unless (defined($overlap_boundary{$x->[$i]->[FID]}))
112            {
113                $overlap_boundary{$x->[$i]->[FID]} = $x->[$i]->[START];
114            }
115    
116            for ($j=$i+1; (($j < @$x) && ($overlap = &entries_overlap($x, $i, $j))); ++$j)
117            {
118                next unless (&same_strand_overlap($x, $i, $j) || &divergent($x, $i, $j));
119    
120                if ($x->[$j]->[STRAND] eq '+')
121                {
122                    $overlap_boundary{$x->[$j]->[FID]}
123                       = &FIG::max( $overlap, $overlap_boundary{$x->[$j]->[FID]} );
124                }
125    
126                if ($x->[$i]->[STRAND] eq '-')
127                {
128                    $overlap_boundary{$x->[$i]->[FID]} = &FIG::max( $overlap, $overlap_boundary{$x->[$i]->[FID]} );
129                }
130            }
131        }
132    }
133    
134    
135    
136    foreach $peg (@pegs)
137    {
138        my $tmp = $peg;
139        if ($peg =~ /^fig\|\d+\.\d+\.peg\.\d+$/)
140        {
141      $genome = $fig->genome_of($peg);      $genome = $fig->genome_of($peg);
142      $loc = $fig->feature_location($peg);      $loc = $fig->feature_location($peg);
143      # print "$len{$peg}\t$peg\t$loc\n";      # print "$len{$peg}\t$peg\t$loc\n";
# Line 25  Line 146 
146          print STDERR "skipping $peg - do not handle complex locs\n";          print STDERR "skipping $peg - do not handle complex locs\n";
147          next;          next;
148      }      }
149        }
150        elsif ($peg =~ /^(new\|(\d+\.\d+)\.peg\.\d+)\t(\S+)/)
151        {
152            $peg = $1;
153            $genome = $2;
154            $loc = $3;
155            ($peg && $genome && $loc) || confess "could not parse \'$tmp\'";
156        }
157        else {
158            confess "Could not parse \'$peg\'";
159        }
160    
161        my $lnC = $fig->contig_ln($genome,$contig);
162      ($contig,$begin,$end) = $fig->boundaries_of($loc);      ($contig,$begin,$end) = $fig->boundaries_of($loc);
163      if ($begin < $end)      if ($begin < $end)
164      {      {
165          $l = sprintf("%s_%d_%d",$contig,1,$begin-1);          $l = sprintf("%s_%d_%d",$contig,1,$begin-1);
166            if ($end < ($lnC-3))
167            {
168          $end += 3;          $end += 3;
169      }      }
170        }
171      else      else
172      {      {
173          $l = sprintf("%s_%d_%d",$contig,$begin+1,$fig->contig_ln($genome,$contig));          $l = sprintf("%s_%d_%d",$contig,$fig->contig_ln($genome,$contig),$begin+1);
174            if ($end > 3)
175            {
176          $end -= 3;          $end -= 3;
177      }      }
178        }
179    
180      $loc = join("_",($contig,$begin,$end));      $loc = join("_",($contig,$begin,$end));
181      $seq = uc $fig->dna_seq($genome,$loc);      $seq = uc $fig->dna_seq($genome,$loc);
182    
183      if (! $seq)      if (! $seq)
184      {      {
185          print STDERR &Dumper($peg,$loc,$seq); die "could not get sequence for $peg";          confess "could not get sequence for fid=$peg, loc=$loc;\n   seq=\'$seq\'";
186      }      }
187    
188  #   The following hack handles uncertainty about whether SEED gives back the STOP or not  #   The following hack handles uncertainty about whether SEED gives back the STOP or not
189      if (substr($seq,-6,3) =~ /TAA|TAG|TGA/i)  
190        if ($fig->translate(substr($seq,-6,3), $code_map) eq qq(*))
191      {      {
192          $end = ($begin < $end) ? $end-3 : $end+3;          $end = ($begin < $end) ? $end-3 : $end+3;
193          $seq = substr($seq,0,length($seq) - 3);          $seq = substr($seq,0,length($seq) - 3);
194      }      }
195      elsif (substr($seq,-3,3) !~ /TAA|TAG|TGA/i)      elsif (($fig->translate(substr($seq,-3,3), $code_map) ne qq(*)) && (not $fig->possibly_truncated($genome, $loc)))
196      {      {
197          die "BAD STOP: $peg $genome $loc $seq\n";          warn "BAD STOP: $peg $genome $loc $seq\n";
198      }      }
199    
     $loc = join("_",($contig,$begin,$end));  
   
200    
201      # BACK INTO THE STOP FROM START OF PEG      # BACK INTO THE STOP FROM START OF PEG
202      $pre_peg = uc $fig->dna_seq($genome,$l);      $pre_peg = uc $fig->dna_seq($genome,$l);
# Line 67  Line 206 
206      for ($i=$len_pre_peg-3; $i > 0 && ! $found; $i-=3)      for ($i=$len_pre_peg-3; $i > 0 && ! $found; $i-=3)
207      {      {
208          $stop =  substr($pre_peg,$i,3);          $stop =  substr($pre_peg,$i,3);
209          if ($stop eq "TAG"  ||  $stop eq "TAA"  ||  $stop eq "TGA")          if ($fig->translate($stop,  $code_map) eq qq(*))
210          {          {
211              # print "FOUND $stop for $l\n";              # print "FOUND $stop for $l\n";
212              $orf = substr($pre_peg,$i+3) . $seq;              $orf = substr($pre_peg,$i+3) . $seq;
213              if (($i-27) > 0)              if (($i-27) > 0)
214              {              {
215                  $pre_orf = substr($pre_peg,$i-27,30);                  $pre_orf = substr($pre_peg,$i-27,30);
216                  my $gs = $fig->org_of($peg);                  my $gs = $fig->genus_species($genome);
217                  my @aliases = grep { $_ =~ /^([NXYZA]P_[0-9\.]+)|(uni\|)|(sp\|)/ } $fig->feature_aliases($peg);                  my @aliases = ();
218                    if ($peg =~ /^fig/)
219                    {
220                        @aliases = grep { $_ =~ /^([NXYZA]P_[0-9\.]+)|(uni\|)|(sp\|)/ } $fig->feature_aliases($peg);
221                    }
222    
223                  print "ID=$peg\n";                  print "ID=$peg\n";
224                  print "GENOME=$gs\n";                  print "GENOME=$gs\n";
225                  if (@aliases > 0)                  if (@aliases > 0)
# Line 83  Line 227 
227                      print "ALIASES=",join(",",@aliases),"\n";                      print "ALIASES=",join(",",@aliases),"\n";
228                  }                  }
229                  print "CONTIG=$contig\n";                  print "CONTIG=$contig\n";
230                  #####print "BEGIN=$begin\n";  
231                  if ($begin < $end)                  if ($begin < $end)
232                  {                  {
233                      print "BEGIN=", $begin-($len_pre_peg-$i-3), "\n";                      print "BEGIN=", $begin-($len_pre_peg-$i-3), "\n";
# Line 96  Line 240 
240                  print "SEQ=$orf\n";                  print "SEQ=$orf\n";
241                  print "PREFIX=$pre_orf\n";                  print "PREFIX=$pre_orf\n";
242                  print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";                  print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";
243                    print "OVERLAP=$overlap_boundary{$peg}\n";
244                  print "///\n";                  print "///\n";
245                  $found = 1;                  $found = 1;
246              }              }
# Line 111  Line 256 
256          print STDERR "DID NOT FIND STOP BEFORE $peg  $contig $begin $end $l\n";          print STDERR "DID NOT FIND STOP BEFORE $peg  $contig $begin $end $l\n";
257      }      }
258  }  }
259    
260    
261    
262    sub from_locus
263    {
264        my ($locus) = @_;
265        my ($contig, $beg, $end) = $fig->boundaries_of($locus);
266        my ($left, $right, $len, $strand);
267    
268        if ($contig && $beg && $end)
269        {
270            $left   = &FIG::min($beg, $end);
271            $right  = &FIG::max($beg, $end);
272            $len    = (1+abs($end-$beg));
273            $strand = ($beg < $end) ? '+' : '-';
274    
275            return ($contig, $beg, $end, $left, $right, $len, $strand);
276        }
277        else
278        {
279            die "Invalid locus $locus";
280        }
281    
282        return ();
283    }
284    
285    sub divergent
286    {
287        my ($x, $i, $j) = @_;
288    
289        my $beg_i   = $x->[$i]->[START];
290        my $end_i   = $x->[$i]->[STOP];
291    
292        my $beg_j   = $x->[$j]->[START];
293        my $end_j   = $x->[$j]->[STOP];
294    
295        if (   &FIG::between($beg_i, $beg_j, $end_i)
296           &&  &FIG::between($beg_j, $beg_i, $end_j)
297           && !&FIG::between($beg_j, $end_i, $end_j)
298           && !&FIG::between($beg_i, $end_j, $end_i)
299           )
300        {
301            return &overlaps($beg_i, $end_i, $beg_j, $end_j);
302        }
303    
304        return 0;
305    }
306    
307    sub same_strand_overlap
308    {
309        my ($x, $i, $j) = @_;
310    
311        my $beg_i   = $x->[$i]->[START];
312        my $end_i   = $x->[$i]->[STOP];
313    
314        my $beg_j   = $x->[$j]->[START];
315        my $end_j   = $x->[$j]->[STOP];
316    
317        if    (  ($x->[$i]->[STRAND] eq '+') && ($x->[$j]->[STRAND] eq '+')
318              && ($beg_i < $beg_j) && ($beg_j <= $end_i) && ($end_i < $end_j)
319              )
320        {
321            return &overlaps($beg_i, $end_i,  $beg_j, $end_j);
322        }
323        elsif (  ($x->[$i]->[STRAND] eq '-') && ($x->[$j]->[STRAND] eq '-')
324              && ($end_i < $end_j) && ($end_j <= $beg_i) && ($beg_i < $beg_j)
325              )
326        {
327            return &overlaps($beg_i, $end_i,  $beg_j, $end_j);
328        }
329        else
330        {
331            &impossible_overlap_warning($x, $i, $j, "Opposite strands in a normal_overlap");
332            return 0;
333        }
334    }
335    
336    sub entries_overlap
337    {
338        my ($x, $i, $j) = @_;
339    
340        return 0 if ($x->[$i]->[CONTIG] ne $x->[$j]->[CONTIG]);
341    
342        $ov = &overlaps( $x->[$i]->[START], $x->[$i]->[STOP]
343                       , $x->[$j]->[START], $x->[$j]->[STOP]);
344    
345        return $ov;
346    }
347    
348    sub overlaps
349    {
350        my ($beg1, $end1, $beg2, $end2) = @_;
351    
352        my ($left1, $right1, $left2, $right2);
353    
354        $left1  = &FIG::min($beg1, $end1);
355        $left2  = &FIG::min($beg2, $end2);
356    
357        $right1 = &FIG::max($beg1, $end1);
358        $right2 = &FIG::max($beg2, $end2);
359    
360        if ($left1 > $left2)
361        {
362            ($left1, $left2)   = ($left2, $left1);
363            ($right1, $right2) = ($right2, $right1);
364        }
365    
366        $ov = 0;
367        if ($right1 >= $left2) { $ov = &FIG::min($right1,$right2) - $left2 + 1; }
368    
369        return $ov;
370    }
371    
372    sub impossible_overlap_warning
373    {
374        my ($x, $i, $j, $msg) = @_;
375    
376        return 0 if $deleted{"$x->[$i]->[FILE]\t$x->[$i]->[LOCUS]"};
377        return 0 if $deleted{"$x->[$j]->[FILE]\t$x->[$j]->[LOCUS]"};
378    
379        if (defined($msg))
380        {
381            print STDERR ("Impossible overlap in $msg:\n"
382                          , join("\t", @ { $x->[$i] }), "\n"
383                          , join("\t", @ { $x->[$j] }), "\n\n")
384                if $ENV{VERBOSE};
385            # confess "$msg";
386        }
387        else
388        {
389            print STDERR ("Impossible overlap:\n$i: "
390                          , join("\t", @ { $x->[$i] }), "\n$j: "
391                          , join("\t", @ { $x->[$j] }), "\n\n")
392                if $ENV{VERBOSE};
393            # confess "aborted";
394        }
395    
396        return 0;
397    }

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.12

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3