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

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.10

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3