[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.12 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3