[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.9, Tue Nov 22 14:57:50 2005 UTC
# Line 1  Line 1 
1  use FIG_Config;  # -*- 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  for ($i=0; $i < @pegs; $i++)  # strip off all but peg ids at front  
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($overlap_boundary{$x->[$i]->[FID]}))
86            {
87                $overlap_boundary{$x->[$i]->[FID]} = $x->[$i]->[START];
88            }
89    
90            for ($j=$i+1; (($j < @$x) && ($overlap = &entries_overlap($x, $i, $j))); ++$j)
91            {
92                next unless (&same_strand_overlap($x, $i, $j) || &divergent($x, $i, $j));
93    
94                if ($x->[$j]->[STRAND] eq '+')
95  {  {
96      if ($pegs[$i] =~ /(\S+)\t.*/)                  $overlap_boundary{$x->[$j]->[FID]}
97                       = &FIG::max( $overlap, $overlap_boundary{$x->[$j]->[FID]} );
98                }
99    
100                if ($x->[$i]->[STRAND] eq '-')
101      {      {
102          $pegs[$i] = $1;                  $overlap_boundary{$x->[$i]->[FID]} = &FIG::max( $overlap, $overlap_boundary{$x->[$i]->[FID]} );
103                }
104      }      }
105  }  }
106    }
107    
108    
109    
110  foreach $peg (@pegs)  foreach $peg (@pegs)
111  {  {
112        if ($peg =~ /^fig\|\d+\.\d+\.peg\.\d+$/)
113        {
114      $genome = $fig->genome_of($peg);      $genome = $fig->genome_of($peg);
115      $loc = $fig->feature_location($peg);      $loc = $fig->feature_location($peg);
116      # print "$len{$peg}\t$peg\t$loc\n";      # print "$len{$peg}\t$peg\t$loc\n";
# Line 26  Line 119 
119          print STDERR "skipping $peg - do not handle complex locs\n";          print STDERR "skipping $peg - do not handle complex locs\n";
120          next;          next;
121      }      }
122        }
123        elsif ($peg =~ /^(new\|(\d+\.\d+)\.peg\.\d+)\t(\S+)/)
124        {
125            $peg = $1;
126            $genome = $2;
127            $loc = $3;
128        }
129    
130        ($genome && $loc) || confess $peg;
131    
132        my $lnC = $fig->contig_ln($genome,$contig);
133      ($contig,$begin,$end) = $fig->boundaries_of($loc);      ($contig,$begin,$end) = $fig->boundaries_of($loc);
134      if ($begin < $end)      if ($begin < $end)
135      {      {
136          $l = sprintf("%s_%d_%d",$contig,1,$begin-1);          $l = sprintf("%s_%d_%d",$contig,1,$begin-1);
137            if ($end < ($lnC-3))
138            {
139          $end += 3;          $end += 3;
140      }      }
141        }
142      else      else
143      {      {
144          $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);
145            if ($end > 3)
146            {
147          $end -= 3;          $end -= 3;
148      }      }
149        }
150    
151      $loc = join("_",($contig,$begin,$end));      $loc = join("_",($contig,$begin,$end));
152      $seq = uc $fig->dna_seq($genome,$loc);      $seq = uc $fig->dna_seq($genome,$loc);
153    
154        if (! $seq)
155        {
156            print STDERR &Dumper($peg,$loc,$seq); die "could not get sequence for $peg";
157        }
158    
159  #   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
160      if (substr($seq,-6,3) =~ /TAA|TAG|TGA/i)      if (substr($seq,-6,3) =~ /TAA|TAG|TGA/i)
161      {      {
# Line 49  Line 164 
164      }      }
165      elsif (substr($seq,-3,3) !~ /TAA|TAG|TGA/i)      elsif (substr($seq,-3,3) !~ /TAA|TAG|TGA/i)
166      {      {
167          die "BAD STOP: $genome $loc $seq\n";          warn "BAD STOP: $peg $genome $loc $seq\n";
168      }      }
169    
170      $loc = join("_",($contig,$begin,$end));      $loc = join("_",($contig,$begin,$end));
# Line 70  Line 185 
185              if (($i-27) > 0)              if (($i-27) > 0)
186              {              {
187                  $pre_orf = substr($pre_peg,$i-27,30);                  $pre_orf = substr($pre_peg,$i-27,30);
188                  my $gs = $fig->org_of($peg);                  my $gs = $fig->genus_species($genome);
189                  my @aliases = grep { $_ =~ /^([NXYZA]P_[0-9\.]+)|(uni\|)|(sp\|)/ } $fig->feature_aliases($peg);                  my @aliases = ();
190                    if ($peg =~ /^fig/)
191                    {
192                        @aliases = grep { $_ =~ /^([NXYZA]P_[0-9\.]+)|(uni\|)|(sp\|)/ } $fig->feature_aliases($peg);
193                    }
194    
195                  print "ID=$peg\n";                  print "ID=$peg\n";
196                  print "GENOME=$gs\n";                  print "GENOME=$gs\n";
197                  if (@aliases > 0)                  if (@aliases > 0)
# Line 79  Line 199 
199                      print "ALIASES=",join(",",@aliases),"\n";                      print "ALIASES=",join(",",@aliases),"\n";
200                  }                  }
201                  print "CONTIG=$contig\n";                  print "CONTIG=$contig\n";
202                  #####print "BEGIN=$begin\n";  
203                  if ($begin < $end)                  if ($begin < $end)
204                  {                  {
205                      print "BEGIN=", $begin-($len_pre_peg-$i-3), "\n";                      print "BEGIN=", $begin-($len_pre_peg-$i-3), "\n";
# Line 92  Line 212 
212                  print "SEQ=$orf\n";                  print "SEQ=$orf\n";
213                  print "PREFIX=$pre_orf\n";                  print "PREFIX=$pre_orf\n";
214                  print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";                  print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";
215                    print "OVERLAP=$overlap_boundary{$peg}\n";
216                  print "///\n";                  print "///\n";
217                  $found = 1;                  $found = 1;
218              }              }
# Line 107  Line 228 
228          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";
229      }      }
230  }  }
231    
232    
233    
234    sub from_locus
235    {
236        my ($locus) = @_;
237        my ($contig, $beg, $end) = $fig->boundaries_of($locus);
238        my ($left, $right, $len, $strand);
239    
240        if ($contig && $beg && $end)
241        {
242            $left   = &FIG::min($beg, $end);
243            $right  = &FIG::max($beg, $end);
244            $len    = (1+abs($end-$beg));
245            $strand = ($beg < $end) ? '+' : '-';
246    
247            return ($contig, $beg, $end, $left, $right, $len, $strand);
248        }
249        else
250        {
251            die "Invalid locus $locus";
252        }
253    
254        return ();
255    }
256    
257    sub divergent
258    {
259        my ($x, $i, $j) = @_;
260    
261        my $beg_i   = $x->[$i]->[START];
262        my $end_i   = $x->[$i]->[STOP];
263    
264        my $beg_j   = $x->[$j]->[START];
265        my $end_j   = $x->[$j]->[STOP];
266    
267        if (   &FIG::between($beg_i, $beg_j, $end_i)
268           &&  &FIG::between($beg_j, $beg_i, $end_j)
269           && !&FIG::between($beg_j, $end_i, $end_j)
270           && !&FIG::between($beg_i, $end_j, $end_i)
271           )
272        {
273            return &overlaps($beg_i, $end_i, $beg_j, $end_j);
274        }
275    
276        return 0;
277    }
278    
279    sub same_strand_overlap
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    (  ($x->[$i]->[STRAND] eq '+') && ($x->[$j]->[STRAND] eq '+')
290              && ($beg_i < $beg_j) && ($beg_j <= $end_i) && ($end_i < $end_j)
291              )
292        {
293            return &overlaps($beg_i, $end_i,  $beg_j, $end_j);
294        }
295        elsif (  ($x->[$i]->[STRAND] eq '-') && ($x->[$j]->[STRAND] eq '-')
296              && ($end_i < $end_j) && ($end_j <= $beg_i) && ($beg_i < $beg_j)
297              )
298        {
299            return &overlaps($beg_i, $end_i,  $beg_j, $end_j);
300        }
301        else
302        {
303            &impossible_overlap_warning($x, $i, $j, "Opposite strands in a normal_overlap");
304            return 0;
305        }
306    }
307    
308    sub entries_overlap
309    {
310        my ($x, $i, $j) = @_;
311    
312        return 0 if ($x->[$i]->[CONTIG] ne $x->[$j]->[CONTIG]);
313    
314        $ov = &overlaps( $x->[$i]->[START], $x->[$i]->[STOP]
315                       , $x->[$j]->[START], $x->[$j]->[STOP]);
316    
317        return $ov;
318    }
319    
320    sub overlaps
321    {
322        my ($beg1, $end1, $beg2, $end2) = @_;
323    
324        my ($left1, $right1, $left2, $right2);
325    
326        $left1  = &FIG::min($beg1, $end1);
327        $left2  = &FIG::min($beg2, $end2);
328    
329        $right1 = &FIG::max($beg1, $end1);
330        $right2 = &FIG::max($beg2, $end2);
331    
332        if ($left1 > $left2)
333        {
334            ($left1, $left2)   = ($left2, $left1);
335            ($right1, $right2) = ($right2, $right1);
336        }
337    
338        $ov = 0;
339        if ($right1 >= $left2) { $ov = &FIG::min($right1,$right2) - $left2 + 1; }
340    
341        return $ov;
342    }
343    
344    sub impossible_overlap_warning
345    {
346        my ($x, $i, $j, $msg) = @_;
347    
348        return 0 if $deleted{"$x->[$i]->[FILE]\t$x->[$i]->[LOCUS]"};
349        return 0 if $deleted{"$x->[$j]->[FILE]\t$x->[$j]->[LOCUS]"};
350    
351        if (defined($msg))
352        {
353            print STDERR ("Impossible overlap in $msg:\n"
354                          , join("\t", @ { $x->[$i] }), "\n"
355                          , join("\t", @ { $x->[$j] }), "\n\n")
356                if $ENV{VERBOSE};
357            # confess "$msg";
358        }
359        else
360        {
361            print STDERR ("Impossible overlap:\n$i: "
362                          , join("\t", @ { $x->[$i] }), "\n$j: "
363                          , join("\t", @ { $x->[$j] }), "\n\n")
364                if $ENV{VERBOSE};
365            # confess "aborted";
366        }
367    
368        return 0;
369    }

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3