[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.7, Sat Aug 6 14:58:56 2005 UTC revision 1.8, Sat Nov 5 22:59:44 2005 UTC
# Line 1  Line 1 
1    # -*- perl -*-
2    
3  use FIG;  use FIG;
4  my $fig = new FIG;  my $fig = new FIG;
5    
6  $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";
7    
8    use constant FID    =>  0;
9    use constant LOCUS  =>  1;
10    use constant CONTIG =>  2;
11    use constant START  =>  3;
12    use constant STOP   =>  4;
13    use constant LEFT   =>  5;
14    use constant RIGHT  =>  6;
15    use constant LEN    =>  7;
16    use constant STRAND =>  8;
17    
18    
19    $tbl_file = "";
20    $trouble  = 0;
21    while (@ARGV)
22    {
23        if ($ARGV[0] =~ m/-h(elp)?/)
24        {
25            die "\n\tusage: $usage\n\n";
26        }
27        elsif (-s $ARGV[0])
28        {
29            $tbl_file = $ARGV[0];
30        }
31        else
32        {
33            warn "Invalid argument $ARGV[0]\n";
34        }
35        shift @ARGV;
36    }
37    die "\nusage: $usage\n\n" if $trouble;
38    
39  @pegs = <STDIN>;  @pegs = <STDIN>;
40  chomp @pegs;  chomp @pegs;
41  # print "NPEGS=", scalar(@pegs), "\n";  # print "NPEGS=", scalar(@pegs), "\n";
42    
43    if ($tbl_file)
44    {
45        open(TBL, "<$tbl_file") || die "Could not read-open $tbl_file";
46        while (defined($entry = <TBL>))
47        {
48            if ($entry =~ m/^(\S+)\t(\S+)/)
49            {
50                ($fid, $loc) = ($1, $2);
51                ($contig, $beg, $end, $left, $right, $len, $strand) = &from_locus($loc);
52                $tbl->{$contig}->{$strand.$end}
53                    = [ $fid, $loc, $contig, $beg, $end, $left, $right, $len, $strand ];
54            }
55            else
56            {
57                print STDERR "Skipping invalid tbl entry: $entry";
58            }
59        }
60    }
61    
62    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
63    #...If there are "new" PEGs, overwrite any existing tbl entries...
64    #-----------------------------------------------------------------------
65    foreach $peg (@pegs)
66    {
67        if ($peg =~ /^(new\|(\d+\.\d+)\.peg\.\d+)\t(\S+)/)
68        {
69            $fid    = $1;
70            $genome = $2;
71            $loc    = $3;
72    
73            ($contig, $beg, $end, $left, $right, $len, $strand) = &from_locus($loc);
74            $tbl->{$contig}->{$strand.$end}
75                = [ $fid, $loc, $contig, $beg, $end, $left, $right, $len, $strand ];
76        }
77    }
78    
79    foreach $contig (keys %$tbl)
80    {
81        @$x = sort { $a->[LEFT] <=> $b->[LEFT] } values % { $tbl->{$contig} };
82    
83        for ($i=0; $i < @$x; ++$i)
84        {
85            unless (defined($overlaps{$x->[$i]->[FID]}))  { $overlaps{$x->[$i]->[FID]} = 0; }
86    
87            for ($j=$i+1; (($j < @$x) && ($overlap = &entries_overlap($x, $i, $j))); ++$j)
88            {
89                next unless (&same_strand_overlap($x, $i, $j) || &divergent($x, $i, $j));
90    
91                if ($x->[$j]->[STRAND] eq '+')
92                {
93                    $overlaps{$x->[$j]->[FID]} = &FIG::max( $overlap, $overlaps{$x->[$j]->[FID]} );
94                }
95    
96                if ($x->[$i]->[STRAND] eq '-')
97                {
98                    $overlaps{$x->[$i]->[FID]} = &FIG::max( $overlap, $overlaps{$x->[$i]->[FID]} );
99                }
100            }
101        }
102    }
103    
104    
105    
106  foreach $peg (@pegs)  foreach $peg (@pegs)
107  {  {
108      if ($peg =~ /^fig\|\d+\.\d+\.peg\.\d+$/)      if ($peg =~ /^fig\|\d+\.\d+\.peg\.\d+$/)
# Line 113  Line 208 
208                  print "SEQ=$orf\n";                  print "SEQ=$orf\n";
209                  print "PREFIX=$pre_orf\n";                  print "PREFIX=$pre_orf\n";
210                  print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";                  print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";
211                    print "OVERLAP=$overlaps{$peg}\n";
212                  print "///\n";                  print "///\n";
213                  $found = 1;                  $found = 1;
214              }              }
# Line 128  Line 224 
224          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";
225      }      }
226  }  }
227    
228    
229    
230    sub from_locus
231    {
232        my ($locus) = @_;
233        my ($contig, $beg, $end) = $fig->boundaries_of($locus);
234        my ($left, $right, $len, $strand);
235    
236        if ($contig && $beg && $end)
237        {
238            $left   = &FIG::min($beg, $end);
239            $right  = &FIG::max($beg, $end);
240            $len    = (1+abs($end-$beg));
241            $strand = ($beg < $end) ? '+' : '-';
242    
243            return ($contig, $beg, $end, $left, $right, $len, $strand);
244        }
245        else
246        {
247            die "Invalid locus $locus";
248        }
249    
250        return ();
251    }
252    
253    sub divergent
254    {
255        my ($x, $i, $j) = @_;
256    
257        my $beg_i   = $x->[$i]->[START];
258        my $end_i   = $x->[$i]->[STOP];
259    
260        my $beg_j   = $x->[$j]->[START];
261        my $end_j   = $x->[$j]->[STOP];
262    
263        if (   &FIG::between($beg_i, $beg_j, $end_i)
264           &&  &FIG::between($beg_j, $beg_i, $end_j)
265           && !&FIG::between($beg_j, $end_i, $end_j)
266           && !&FIG::between($beg_i, $end_j, $end_i)
267           )
268        {
269            return &overlaps($beg_i, $end_i, $beg_j, $end_j);
270        }
271    
272        return 0;
273    }
274    
275    sub same_strand_overlap
276    {
277        my ($x, $i, $j) = @_;
278    
279        my $beg_i   = $x->[$i]->[START];
280        my $end_i   = $x->[$i]->[STOP];
281    
282        my $beg_j   = $x->[$j]->[START];
283        my $end_j   = $x->[$j]->[STOP];
284    
285        if    (  ($x->[$i]->[STRAND] eq '+') && ($x->[$j]->[STRAND] eq '+')
286              && ($beg_i < $beg_j) && ($beg_j <= $end_i) && ($end_i < $end_j)
287              )
288        {
289            return &overlaps($beg_i, $end_i,  $beg_j, $end_j);
290        }
291        elsif (  ($x->[$i]->[STRAND] eq '-') && ($x->[$j]->[STRAND] eq '-')
292              && ($end_i < $end_j) && ($end_j <= $beg_i) && ($beg_i < $beg_j)
293              )
294        {
295            return &overlaps($beg_i, $end_i,  $beg_j, $end_j);
296        }
297        else
298        {
299            &impossible_overlap_warning($x, $i, $j, "Opposite strands in a normal_overlap");
300            return 0;
301        }
302    }
303    
304    sub entries_overlap
305    {
306        my ($x, $i, $j) = @_;
307    
308        return 0 if ($x->[$i]->[CONTIG] ne $x->[$j]->[CONTIG]);
309    
310        $ov = &overlaps( $x->[$i]->[START], $x->[$i]->[STOP]
311                       , $x->[$j]->[START], $x->[$j]->[STOP]);
312    
313        return $ov;
314    }
315    
316    sub overlaps
317    {
318        my ($beg1, $end1, $beg2, $end2) = @_;
319    
320        my ($left1, $right1, $left2, $right2);
321    
322        $left1  = &FIG::min($beg1, $end1);
323        $left2  = &FIG::min($beg2, $end2);
324    
325        $right1 = &FIG::max($beg1, $end1);
326        $right2 = &FIG::max($beg2, $end2);
327    
328        if ($left1 > $left2)
329        {
330            ($left1, $left2)   = ($left2, $left1);
331            ($right1, $right2) = ($right2, $right1);
332        }
333    
334        $ov = 0;
335        if ($right1 >= $left2) { $ov = &FIG::min($right1,$right2) - $left2 + 1; }
336    
337        return $ov;
338    }
339    
340    sub impossible_overlap_warning
341    {
342        my ($x, $i, $j, $msg) = @_;
343    
344        return 0 if $deleted{"$x->[$i]->[FILE]\t$x->[$i]->[LOCUS]"};
345        return 0 if $deleted{"$x->[$j]->[FILE]\t$x->[$j]->[LOCUS]"};
346    
347        if (defined($msg))
348        {
349            print STDERR ("Impossible overlap in $msg:\n"
350                          , join("\t", @ { $x->[$i] }), "\n"
351                          , join("\t", @ { $x->[$j] }), "\n\n")
352                if $ENV{VERBOSE};
353            # confess "$msg";
354        }
355        else
356        {
357            print STDERR ("Impossible overlap:\n$i: "
358                          , join("\t", @ { $x->[$i] }), "\n$j: "
359                          , join("\t", @ { $x->[$j] }), "\n\n")
360                if $ENV{VERBOSE};
361            # confess "aborted";
362        }
363    
364        return 0;
365    }

Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3