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

Annotation of /FigKernelScripts/build_initial_objects_for_start_predictions.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (view) (download) (as text)

1 : overbeek 1.8 # -*- perl -*-
2 : overbeek 1.1
3 :     use FIG;
4 :     my $fig = new FIG;
5 :    
6 : overbeek 1.8 $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 : overbeek 1.1
39 :     @pegs = <STDIN>;
40 :     chomp @pegs;
41 :     # print "NPEGS=", scalar(@pegs), "\n";
42 : overbeek 1.5
43 : overbeek 1.8 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 : overbeek 1.5 foreach $peg (@pegs)
107 : overbeek 1.1 {
108 : overbeek 1.7 if ($peg =~ /^fig\|\d+\.\d+\.peg\.\d+$/)
109 : overbeek 1.1 {
110 : overbeek 1.5 $genome = $fig->genome_of($peg);
111 :     $loc = $fig->feature_location($peg);
112 :     # print "$len{$peg}\t$peg\t$loc\n";
113 :     if ($loc =~ /,/)
114 :     {
115 :     print STDERR "skipping $peg - do not handle complex locs\n";
116 :     next;
117 :     }
118 : overbeek 1.1 }
119 : overbeek 1.5 elsif ($peg =~ /^(new\|(\d+\.\d+)\.peg\.\d+)\t(\S+)/)
120 :     {
121 :     $peg = $1;
122 :     $genome = $2;
123 :     $loc = $3;
124 : overbeek 1.1 }
125 :    
126 : overbeek 1.5 ($genome && $loc) || confess $peg;
127 :    
128 : overbeek 1.4 my $lnC = $fig->contig_ln($genome,$contig);
129 : overbeek 1.1 ($contig,$begin,$end) = $fig->boundaries_of($loc);
130 :     if ($begin < $end)
131 :     {
132 :     $l = sprintf("%s_%d_%d",$contig,1,$begin-1);
133 : overbeek 1.4 if ($end < ($lnC-3))
134 :     {
135 :     $end += 3;
136 :     }
137 : overbeek 1.1 }
138 :     else
139 :     {
140 : overbeek 1.6 $l = sprintf("%s_%d_%d",$contig,$fig->contig_ln($genome,$contig),$begin+1);
141 : overbeek 1.4 if ($end > 3)
142 :     {
143 :     $end -= 3;
144 :     }
145 : overbeek 1.1 }
146 :    
147 :     $loc = join("_",($contig,$begin,$end));
148 :     $seq = uc $fig->dna_seq($genome,$loc);
149 : overbeek 1.5
150 : overbeek 1.3 if (! $seq)
151 :     {
152 :     print STDERR &Dumper($peg,$loc,$seq); die "could not get sequence for $peg";
153 :     }
154 :    
155 : overbeek 1.1 # The following hack handles uncertainty about whether SEED gives back the STOP or not
156 :     if (substr($seq,-6,3) =~ /TAA|TAG|TGA/i)
157 :     {
158 : overbeek 1.2 $end = ($begin < $end) ? $end-3 : $end+3;
159 : overbeek 1.1 $seq = substr($seq,0,length($seq) - 3);
160 :     }
161 :     elsif (substr($seq,-3,3) !~ /TAA|TAG|TGA/i)
162 :     {
163 : overbeek 1.4 warn "BAD STOP: $peg $genome $loc $seq\n";
164 : overbeek 1.1 }
165 :    
166 : overbeek 1.2 $loc = join("_",($contig,$begin,$end));
167 : overbeek 1.1
168 :    
169 :     # BACK INTO THE STOP FROM START OF PEG
170 :     $pre_peg = uc $fig->dna_seq($genome,$l);
171 :     $found = 0;
172 :     $len_pre_peg = length($pre_peg);
173 :    
174 :     for ($i=$len_pre_peg-3; $i > 0 && ! $found; $i-=3)
175 :     {
176 :     $stop = substr($pre_peg,$i,3);
177 :     if ($stop eq "TAG" || $stop eq "TAA" || $stop eq "TGA")
178 :     {
179 :     # print "FOUND $stop for $l\n";
180 :     $orf = substr($pre_peg,$i+3) . $seq;
181 :     if (($i-27) > 0)
182 :     {
183 :     $pre_orf = substr($pre_peg,$i-27,30);
184 : overbeek 1.5 my $gs = $fig->genus_species($genome);
185 :     my @aliases = ();
186 :     if ($peg =~ /^fig/)
187 :     {
188 :     @aliases = grep { $_ =~ /^([NXYZA]P_[0-9\.]+)|(uni\|)|(sp\|)/ } $fig->feature_aliases($peg);
189 :     }
190 :    
191 : overbeek 1.1 print "ID=$peg\n";
192 :     print "GENOME=$gs\n";
193 :     if (@aliases > 0)
194 :     {
195 :     print "ALIASES=",join(",",@aliases),"\n";
196 :     }
197 :     print "CONTIG=$contig\n";
198 : overbeek 1.5
199 : overbeek 1.1 if ($begin < $end)
200 :     {
201 :     print "BEGIN=", $begin-($len_pre_peg-$i-3), "\n";
202 :     }
203 :     else
204 :     {
205 :     print "BEGIN=", $begin+($len_pre_peg-$i-3), "\n";
206 :     }
207 :     print "END=$end\n";
208 :     print "SEQ=$orf\n";
209 :     print "PREFIX=$pre_orf\n";
210 :     print "OLD_START_POS=", $len_pre_peg-$i-2, "\n";
211 : overbeek 1.8 print "OVERLAP=$overlaps{$peg}\n";
212 : overbeek 1.1 print "///\n";
213 :     $found = 1;
214 :     }
215 :     else
216 :     {
217 :     print STDERR "skipping $peg - not enough prefix before orf\n";
218 :     }
219 :     last;
220 :     }
221 :     }
222 :     if ( ! $found)
223 :     {
224 :     print STDERR "DID NOT FIND STOP BEFORE $peg $contig $begin $end $l\n";
225 :     }
226 :     }
227 : overbeek 1.8
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 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3