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

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.11

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3