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

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.14 ########################################################################
2 : olson 1.12 #
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 : efrank 1.1 use FIG;
20 :    
21 : efrank 1.2 $usage = "usage: parse_genbank Genome Dir < genbank.entry";
22 : gdpusch 1.9 while ($ARGV[0] =~ m/^-/)
23 :     {
24 :     if ($ARGV[0] =~ m/^-{1,2}genome=\S+/)
25 :     {
26 :     $force_genome = shift @ARGV;
27 :     $force_genome =~ s/^-{1,2}genome=\"?//;
28 :     $force_genome =~ s/\"$//;
29 :     $force_genome =~ s/\s+/ /g;
30 :    
31 :     print STDERR "\nGenome bioname will be taken as: \"$force_genome\"\n";
32 :     }
33 :     elsif ($ARGV[0] =~ m/^-{1,2}tax=\S+/)
34 :     {
35 :     $force_tax = shift @ARGV;
36 :     $force_tax =~ s/^-{1,2}tax=\"?//;
37 :     $force_tax =~ s/\"$//;
38 :     $force_tax =~ s/\s+/ /g;
39 :     if ($force_tax !~ m/\.$/) { $force_tax .= "."; }
40 :    
41 :     print STDERR "\nTaxonomy will be taken as: \"$force_tax\"\n";
42 :     }
43 :     else
44 :     {
45 :     die "Could not handle $ARGV[0]";
46 :     }
47 :     }
48 :    
49 : efrank 1.2 (
50 : efrank 1.3 ($genome = shift(@ARGV)) &&
51 :     ($dir = shift(@ARGV))
52 : efrank 1.1 )
53 :     || die $usage;
54 :    
55 : efrank 1.2 $prefixP = "fig|$genome.peg.";
56 :     $prefixR = "fig|$genome.rna.";
57 :    
58 : efrank 1.3 &FIG::verify_dir("$dir/Features/peg");
59 :     &FIG::verify_dir("$dir/Features/rna");
60 :     open(CONTIGS,">$dir/contigs") || die "could not open $dir/contigs";
61 : efrank 1.2 open(TBLPEG,">$dir/Features/peg/tbl") || die "could not open $dir/Features/peg/tbl";
62 :     open(FASTAPEG,">$dir/Features/peg/fasta") || die "could not open $dir/Features/peg/fasta";
63 : efrank 1.3 open(ASSIGNMENTS,">$dir/assigned_functions") || die "could not open $dir/assigned_functions";
64 : efrank 1.2 open(TBLRNA,">$dir/Features/rna/tbl") || die "could not open $dir/Features/rna/tbl";
65 :     open(FASTARNA,">$dir/Features/rna/fasta") || die "could not open $dir/Features/rna/fasta";
66 :    
67 : overbeek 1.4 open(TAX,">$dir/TAXONOMY") || die "could not open $dir/TAXONOMY";
68 :     open(GENOME,">$dir/GENOME") || die "could not open $dir/GENOME";
69 :     open(PROJECT,">$dir/PROJECT") || die "could not open $dir/PROJECT";
70 :     print PROJECT "Parsed NCBI Reference Sequences\n";
71 :     close(PROJECT);
72 :    
73 :     $idNp = 1;
74 :     $idNr = 1;
75 :    
76 : gdpusch 1.9 $/ = "\n//\n";
77 :     while ($record = <STDIN>)
78 : efrank 1.1 {
79 : gdpusch 1.9 if ($record !~ /LOCUS\s+(\S+)/s)
80 :     {
81 :    
82 :     print STDERR "No LOCUS line for record begining with:\n"
83 :     , substr($record, 0, 160), "\n\n";
84 :     }
85 :     else
86 : efrank 1.1 {
87 :     $id = $1;
88 : gdpusch 1.9 if ($record =~ /ORIGIN(.*?)(\/\/|LOCUS)/s)
89 : efrank 1.1 {
90 :     $seq = $1;
91 :     $seq =~ s/\s//gs;
92 :     $seq =~ s/\d//g;
93 :     undef $contigs;
94 :     $contigs->{$id} = $seq;
95 :     &FIG::display_id_and_seq($id,\$seq,\*CONTIGS);
96 :     }
97 :     else
98 :     {
99 : gdpusch 1.9 warn "could not find contig sequence for $id\n";
100 : overbeek 1.4 next;
101 :     }
102 :    
103 : overbeek 1.14 if ($record =~ /\n {0,4}ORGANISM\s+(\S[^\n]+(\n\s{10,14}\S[^\n]+)*)/s)
104 : overbeek 1.4 {
105 : overbeek 1.14 my $block = $1;
106 :     my @lines = split(/\n/,$block);
107 :    
108 :     my @genome = ();
109 :     my @full_tax = ();
110 :     for ($i=0; ($i < @lines) && ($lines[$i] !~ /;/); $i++)
111 :     {
112 :     push(@genome,$lines[$i]);
113 :     }
114 :     while ($i < @lines)
115 :     {
116 :     push(@full_tax,$lines[$i]);
117 :     $i++;
118 :     }
119 :    
120 :     $genome = join(" ",map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @genome);
121 :     $tax = join(" ",map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @full_tax);
122 :    
123 : overbeek 1.4 $tax =~ s/\n\s+//g;
124 : overbeek 1.8 $tax =~ s/ {2,}/ /g;
125 : overbeek 1.4 if (! $written_genome)
126 :     {
127 : gdpusch 1.9 if ($force_genome) { $genome = $force_genome; }
128 : overbeek 1.4 print GENOME "$genome\n";
129 :     close(GENOME);
130 : gdpusch 1.9
131 :     if ($force_tax) { $tax = $force_tax; }
132 : overbeek 1.4 print TAX "$tax\n";
133 :     close(TAX);
134 :    
135 :     $written_genome = "$tax\t$genome";
136 :     }
137 : gdpusch 1.9 elsif (($written_genome ne "$tax\t$genome") && (!$force_genome || !$force_tax))
138 : overbeek 1.4 {
139 :     print STDERR "serious mismatch in GENOME/TAX for $id\n$written_genome\n$tax\t$genome\n\n";
140 :     }
141 : efrank 1.1 }
142 :    
143 : gdpusch 1.9 while ($record =~ /\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/gs)
144 : efrank 1.1 {
145 :     $cds = $1;
146 : overbeek 1.8 if (($cds !~ /\/pseudo/) && (($cds !~ /\/exception/) || ($cds =~ /\/translation/)))
147 : overbeek 1.7 {
148 :     &process_cds($id,\$cds,$prefixP,\$idNp,$contigs,\*TBLPEG,\*FASTAPEG,\*ASSIGNMENTS);
149 :     }
150 : efrank 1.2 }
151 :    
152 : gdpusch 1.9 while ($record =~ /\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/gs)
153 : efrank 1.2 {
154 :     $rna = $1;
155 :     &process_rna($id,\$rna,$prefixR,\$idNr,$contigs,\*TBLRNA,\*FASTARNA);
156 : efrank 1.1 }
157 :     }
158 :     }
159 : efrank 1.2 close(CONTIGS);
160 :     close(TBLPEG);
161 : overbeek 1.4 close(FASTAPEG);
162 : efrank 1.2 close(ASSIGNMENTS);
163 :     close(TBLRNA);
164 :     close(FASTARNA);
165 :    
166 :     if (! -s "$dir/contigs") { unlink("$dir/contigs"); print STDERR "no contigs in $dir\n"; }
167 : efrank 1.3 if (! -s "$dir/assigned_functions") { unlink("$dir/assigned_functions"); print STDERR "no assigned_functions in $dir\n"; }
168 :     if (! -s "$dir/Features/peg/tbl") { system "rm -rf $dir/Features/peg"; print STDERR "no PEGs in $dir\n"; }
169 :     if (! -s "$dir/Features/rna/tbl") { system "rm -rf $dir/Features/rna"; print STDERR "no RNAs in $dir\n"; }
170 : overbeek 1.8 if ((! -s "$dir/contigs") && (! -s "$dir/Features/peg/tbl")) { system "rm -rf $dir" }
171 : efrank 1.2
172 : efrank 1.1
173 :     sub process_cds {
174 : overbeek 1.4 my($contigID,$cdsP,$prefix,$idNp,$contigs,$fh_tbl,$fh_fasta,$fh_ass) = @_;
175 : gdpusch 1.9 my($loc,@aliases,$func,$trans,$id,$precise,$dna,$prot);
176 :    
177 : overbeek 1.13 ++$recnum;
178 : efrank 1.1 ($loc,$precise) = &get_loc($contigID,$cdsP);
179 :     @aliases = &get_aliases($cdsP);
180 :     $func = &get_func($cdsP);
181 :     $trans = &get_trans($cdsP);
182 :    
183 :     if ($loc || $trans)
184 :     {
185 : gdpusch 1.9 if ($dna = &FIG::extract_seq($contigs,$loc))
186 : efrank 1.1 {
187 : gdpusch 1.9 if ($precise)
188 : efrank 1.1 {
189 : gdpusch 1.9 $prot = &FIG::translate($dna,undef,1);
190 :     if ($trans && (uc $prot ne uc $trans))
191 :     {
192 : overbeek 1.4 # print STDERR "BAD TRANSLATION: $contigID $$cdsP\n";
193 :     # print STDERR &Dumper($trans,$prot,$dna),"\n";
194 : gdpusch 1.9 }
195 : efrank 1.1 }
196 :     }
197 : gdpusch 1.9 else
198 :     {
199 :     warn "Could not get DNA sequence for CDS at $loc for entry:\n$$cdsP\nof record begining with:\n"
200 :     , substr($record, 0, 160), "\n\n";
201 :     }
202 :    
203 : overbeek 1.4 $id = $prefix . "$$idNp";
204 :     $$idNp++;
205 : efrank 1.1 print $fh_tbl "$id\t$loc\t",join("\t",@aliases),"\n";
206 : overbeek 1.4
207 : efrank 1.1 if ($func)
208 :     {
209 :     print $fh_ass "$id\t$func\n";
210 :     }
211 : overbeek 1.13
212 : efrank 1.1 if ($trans)
213 :     {
214 :     &FIG::display_id_and_seq($id,\$trans,$fh_fasta);
215 :     }
216 : gdpusch 1.9 else
217 :     {
218 :     warn "No translation for $id";
219 :     }
220 : efrank 1.1 }
221 :     else
222 :     {
223 :     print &Dumper($cdsP); die "aborted";
224 :     }
225 :     }
226 :    
227 :     sub get_loc {
228 :     my($contigID,$cdsP) = @_;
229 :     my($beg,$end,$loc,$locS,$iter,$n,$args,@pieces);
230 :     my($precise);
231 :    
232 : efrank 1.2 if ($$cdsP =~ /^(([^\n]+)((\n\s{21,23}[^\/ ][^\n]+)+)?)/s)
233 : efrank 1.1 {
234 :     $locS = $1;
235 :     $locS =~ s/\s//g;
236 :     $precise = ($locS !~ /[<>]/);
237 : gdpusch 1.10
238 : efrank 1.1 @pieces = ();
239 :     $n = 0;
240 :     $iter = 500;
241 :     while (($locS !~ /^\[\d+\]$/) && $iter)
242 :     {
243 :     $iter--;
244 :     if ($locS =~ s/[<>]?(\d+)\.\.[<>]?(\d+)/\[$n\]/)
245 :     {
246 :     push(@pieces,["loc",$1,$2]);
247 :     $n++;
248 :     }
249 : gdpusch 1.10
250 : efrank 1.1 if ($locS =~ s/([,\(])(\d+)([,\)])/$1\[$n\]$3/)
251 :     {
252 :     push(@pieces,["loc",$2,$2]);
253 :     $n++;
254 :     }
255 : gdpusch 1.10
256 : efrank 1.1 if ($locS =~ s/join\((\[\d+\](,\[\d+\])*)\)/\[$n\]/)
257 :     {
258 :     $args = $1;
259 :     push(@pieces,["join",map { $_ =~ /^\[(\d+)\]$/; $1 } split(/,/,$args)]);
260 :     $n++;
261 :     }
262 : gdpusch 1.10
263 : efrank 1.1 if ($locS =~ s/complement\((\[\d+\](,\[\d+\])*)\)/\[$n\]/g)
264 :     {
265 :     $args = $1;
266 :     push(@pieces,["complement",map { $_ =~ /^\[(\d+)\]$/; $1 } split(/,/,$args)]);
267 :     $n++;
268 :     }
269 :     }
270 : gdpusch 1.10
271 : efrank 1.1 if ($locS =~ /^\[(\d+)\]$/)
272 :     {
273 :     $loc = &conv(\@pieces,$1,$contigID);
274 :     }
275 :     else
276 :     {
277 :     print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
278 :     die "aborted";
279 :     }
280 : gdpusch 1.10
281 : efrank 1.1 my @locs = split(/,/,$loc);
282 : gdpusch 1.10 #...STOP is now included, so don't trim it off...
283 :     # &trim_stop(\@locs);
284 : efrank 1.1 $loc = join(",",@locs);
285 :     }
286 :     else
287 :     {
288 :     print STDERR &Dumper($cdsP); die "could not parse location";
289 :     die "aborted";
290 :     }
291 :     return ($loc,$precise);
292 :     }
293 :    
294 :     sub trim_stop {
295 :     my($locs) = @_;
296 :     my($to_trim,$n);
297 :    
298 :     $to_trim = 3;
299 :     while ($to_trim && (@$locs > 0))
300 :     {
301 :     $n = @$locs - 1;
302 :     if ($locs->[$n] =~ /^(\S+)_(\d+)_(\d+)$/)
303 :     {
304 :     if ($2 <= ($3-$to_trim))
305 :     {
306 :     $locs->[$n] = join("_",($1,$2,$3-$to_trim));
307 :     $to_trim = 0;
308 :     }
309 :     elsif ($2 >= ($3 + $to_trim))
310 :     {
311 :     $locs->[$n] = join("_",($1,$2,$3+$to_trim));
312 :     $to_trim = 0;
313 :     }
314 :     else
315 :     {
316 :     $to_trim -= abs($3-$2) + 1;
317 :     pop @$locs;
318 :     }
319 :     }
320 :     else
321 :     {
322 :     die "could not parse $locs->[$n]";
323 :     }
324 :     }
325 :     }
326 :    
327 :    
328 :     sub conv {
329 :     my($pieces,$n,$contigID) = @_;
330 :    
331 :     if ($pieces->[$n]->[0] eq "loc")
332 :     {
333 :     return join("_",$contigID,@{$pieces->[$n]}[1..2]);
334 :     }
335 :     elsif ($pieces->[$n]->[0] eq "join")
336 :     {
337 :     return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
338 :     }
339 :     elsif ($pieces->[$n]->[0] eq "complement")
340 :     {
341 :     return join(",",&complement(join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
342 :     }
343 :     }
344 :    
345 :     sub complement {
346 :     my($loc) = @_;
347 :     my(@locs);
348 :    
349 :     @locs = reverse split(/,/,$loc);
350 :     foreach $loc (@locs)
351 :     {
352 :     $loc =~ /^(\S+)_(\d+)_(\d+)$/;
353 :     $loc = join("_",($1,$3,$2));
354 :     }
355 :     return join(",",@locs);
356 :     }
357 :    
358 :     sub get_aliases {
359 :     my($cdsP) = @_;
360 :     my($db_ref);
361 :    
362 :     my @aliases = ();
363 : overbeek 1.4 while ($$cdsP =~ /\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/sg)
364 : efrank 1.1 {
365 :     push(@aliases,$2);
366 :     }
367 :    
368 : overbeek 1.4 while ($$cdsP =~ /\/db_xref=\"([^\"]+)\"/sg)
369 : efrank 1.1 {
370 :     $db_ref = $1;
371 :     $db_ref =~ s/^GI:/gi\|/;
372 : overbeek 1.4 $db_ref =~ s/SWISS-PROT:/sp\|/;
373 : efrank 1.1 push(@aliases,$db_ref);
374 :     }
375 :    
376 :     return @aliases;
377 :     }
378 :    
379 :     sub get_trans {
380 :     my($cdsP) = @_;
381 :     my $tran = "";
382 :    
383 :     if ($$cdsP =~ /\/translation=\"([^"]*)\"/s)
384 :     {
385 :     $tran = $1;
386 :     $tran =~ s/\s//gs;
387 :     }
388 : overbeek 1.8 elsif ($$cdsP =~ /\/protein_id=\"([^"]+)\"/)
389 :     {
390 :     $tran = $1;
391 :     $tran =~ s/\s//gs;
392 :     }
393 : efrank 1.1 return $tran;
394 :     }
395 :    
396 :     sub get_func {
397 :     my($cdsP) = @_;
398 :     my $func = "";
399 : overbeek 1.11
400 : overbeek 1.13 print STDERR "\nRecord $recnum:\n$$cdsP\n" if $ENV{VERBOSE};
401 :     if (($$cdsP =~ /\/function=\"([^"]*)\"/s) && ($func = $1) && &ok_func($func))
402 : overbeek 1.11 {
403 : overbeek 1.13 print STDERR "Branch 1: $func\n" if $ENV{VERBOSE};
404 : overbeek 1.11 }
405 : overbeek 1.13 elsif (($$cdsP =~ /\/product=\"([^"]*)\"/s) && ($func = $1) && &ok_func($func)
406 :     && (($func =~ / /) || ($$cdsP !~ /\/note/)))
407 : overbeek 1.4 {
408 : overbeek 1.13 print STDERR "Branch 2: $func\n" if $ENV{VERBOSE};
409 : overbeek 1.4 }
410 : overbeek 1.13 elsif (($$cdsP =~ /\/note=\"([^"]*)\"/s) && ($func = $1) && &ok_func($func))
411 : overbeek 1.4 {
412 : overbeek 1.13 print STDERR "Branch 3: $func\n" if $ENV{VERBOSE};
413 : overbeek 1.4 }
414 : overbeek 1.11 else
415 : efrank 1.2 {
416 : overbeek 1.13 print STDERR "No non-hypo found\n" if $ENV{VERBOSE};
417 : efrank 1.2 }
418 : overbeek 1.11 $func =~ s/\s+/ /gs;
419 : overbeek 1.13 print STDERR " --> $func\n" if $ENV{VERBOSE};
420 : overbeek 1.11
421 : overbeek 1.4 $func = &fixup_func($func);
422 : overbeek 1.13 print STDERR "Retval: $func\n\n" if $ENV{VERBOSE};
423 : overbeek 1.11
424 : overbeek 1.4 return $func;
425 :     }
426 :    
427 :     sub fixup_func {
428 :     my($func) = @_;
429 :    
430 :     $func =~ s/^COG\d+:\s+//i;
431 :     $func =~ s/^ORFID:\S+\s+//i;
432 : efrank 1.2 return $func;
433 :     }
434 :    
435 : overbeek 1.4 sub ok_func {
436 :     my($func) = @_;
437 :    
438 :     return (
439 :     ($func !~ /^[a-zA-Z]{1,3}\d+[a-zA-Z]{0,3}(\.\d+([a-zA-Z])?)?$/) &&
440 :     ($func !~ /^\d+$/) &&
441 :     ($func !~ /ensangp/i) &&
442 :     ($func !~ /agr_/i) &&
443 :     ($func !~ /^SC.*:/) &&
444 :     ($func !~ /^RIKEN/) &&
445 :     ($func !~ /\slen: \d+/) &&
446 :     ($func !~ /^similar to/i) &&
447 :     ($func !~ /^CPX/i) &&
448 :     ($func !~ /^\S{3,4}( or \S+ protein)?$/i) &&
449 :     ($func !~ /^putative$/) &&
450 :     ($func !~ /\scomes from/i) &&
451 :     ($func !~ /^unknown( function)?$/i) &&
452 :     ($func !~ /^hypothetical( (function|protein))?$/i) &&
453 :     ($func !~ /^orf /i) &&
454 :     ($func !~ /^ORF\s?\d+[lr]?$/i) &&
455 :     ($func !~ /^[a-z]{1,3}\s?\d+$/i)
456 :     )
457 :     }
458 :    
459 : efrank 1.2 sub process_rna {
460 : overbeek 1.4 my($contigID,$cdsP,$prefix,$idNr,$contigs,$fh_tbl,$fh_dna) = @_;
461 : gdpusch 1.9 my($loc,@aliases,$func,$trans,$id,$precise,$DNA);
462 :    
463 : efrank 1.2 ($loc,$precise) = &get_loc($contigID,$cdsP);
464 :     $func = &get_rna_func($cdsP);
465 :    
466 :     if ($loc)
467 :     {
468 : gdpusch 1.9 if ($dna = &FIG::extract_seq($contigs,$loc))
469 : efrank 1.2 {
470 : gdpusch 1.9 $id = $prefix . "$$idNr";
471 :     $$idNr++;
472 :     if (! $func )
473 :     {
474 :     warn "could not get func $$cdsP\n";
475 :     }
476 :     else
477 :     {
478 :     print $fh_tbl "$id\t$loc\t$func\n";
479 :     &FIG::display_id_and_seq($id,\$dna,$fh_dna);
480 :     }
481 : efrank 1.2 }
482 :     else
483 :     {
484 : gdpusch 1.9 warn "Could not get DNA sequence for RNA at $loc for entry\n$$cdsP\nof record begining with:\n"
485 :     , substr($record, 0, 160), "\n\n";
486 : efrank 1.2 }
487 :     }
488 :     else
489 :     {
490 :     print &Dumper($cdsP); die "aborted";
491 :     }
492 :     }
493 :    
494 :     sub get_rna_func {
495 :     my($cdsP) = @_;
496 :     my $func = "";
497 :    
498 :     if ($$cdsP =~ /\/product=\"([^"]*)\"/s)
499 :     {
500 :     $func = $1;
501 :     $func =~ s/\s+/ /gs;
502 :     }
503 :     elsif ($$cdsP =~ /\/gene=\"([^"]*)\"/s)
504 :     {
505 :     $func = $1;
506 :     $func =~ s/\s+/ /gs;
507 :     }
508 :     elsif ($$cdsP =~ /\/note=\"([^"]*)\"/s)
509 : efrank 1.1 {
510 :     $func = $1;
511 :     $func =~ s/\s+/ /gs;
512 :     }
513 :     return $func;
514 :     }
515 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3