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

Diff of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.17, Wed May 14 15:56:47 2008 UTC revision 1.18, Fri Jan 30 13:03:27 2009 UTC
# Line 1  Line 1 
1    # -*- Perl -*-
2    
3  ########################################################################  ########################################################################
4  #  #
5  # Copyright (c) 2003-2006 University of Chicago and Fellowship  # Copyright (c) 2003-2006 University of Chicago and Fellowship
# Line 15  Line 17 
17  # Genomes at veronika@thefig.info or download a copy from  # Genomes at veronika@thefig.info or download a copy from
18  # http://www.theseed.org/LICENSE.TXT.  # http://www.theseed.org/LICENSE.TXT.
19  #  #
20    ########################################################################
21    
22  use FIG;  use FIG;
23  use Carp;  use Carp;
24    use Data::Dumper;
25    
26    my $fig = FIG->new();
27    my %trans_table;
28    
29    my $usage = qq(usage: parse_genbank [-bioname=GENOME_BIONAME] [-taxonomy=TAXONOMY] [-source=DATA_SOURCE] Taxon_ID Org_Dir < genbank.entry);
30    
31  $usage = "usage: parse_genbank Genome Dir < genbank.entry";  my $trouble = 0;
32    my $force_bioname;
33    my $force_taxonomy;
34    my $data_source = qq(Parsed GenBank File\n);
35  while ($ARGV[0] =~ m/^-/)  while ($ARGV[0] =~ m/^-/)
36  {  {
37      if    ($ARGV[0] =~ m/^-{1,2}genome=\S+/)      if    ($ARGV[0] =~ m/^-{1,2}bioname=\S+/) {
38      {          $force_bioname =  shift @ARGV;
39          $force_genome =  shift @ARGV;          $force_bioname =~ s/^-{1,2}bioname=\"?//;
40          $force_genome =~ s/^-{1,2}genome=\"?//;          $force_bioname =~ s/\"$//;
41          $force_genome =~ s/\"$//;          $force_bioname =~ s/\s+/ /g;
         $force_genome =~ s/\s+/ /g;  
42    
43          print STDERR "\nGenome bioname will be taken as: \"$force_genome\"\n";          print STDERR qq(\nGenome bioname will be taken as: \"$force_bioname\"\n);
44      }      }
45      elsif ($ARGV[0] =~ m/^-{1,2}tax=\S+/)      elsif ($ARGV[0] =~ m/^-{1,2}taxonomy=\S+/) {
46      {          $force_taxonomy =  shift @ARGV;
47          $force_tax =  shift @ARGV;          $force_taxonomy =~ s/^-{1,2}taxonomy=\"?//;
48          $force_tax =~ s/^-{1,2}tax=\"?//;          $force_taxonomy =~ s/\"$//;
49          $force_tax =~ s/\"$//;          $force_taxonomy =~ s/\s+/ /g;
50          $force_tax =~ s/\s+/ /g;          if ($force_taxonomy !~ m/\.$/) { $force_taxonomy .= qq(.); }
         if ($force_tax !~ m/\.$/) { $force_tax .= "."; }  
51    
52          print STDERR "\nTaxonomy will be taken as: \"$force_tax\"\n";          print STDERR qq(\nTaxonomy will be taken as: \"$force_taxonomy\"\n);
53      }      }
54      else      elsif ($ARGV[0] =~ m/^-{1,2}source=\S+/) {
55      {          $data_source =  shift @ARGV;
56          die "Could not handle $ARGV[0]";          $data_source =~ s/^-{1,2}source=[\"\']//;
57            $data_source =~ s/[\"\']$//;
58    
59            print STDERR qq(Data source string is: \"$data_source\"\n);
60        }
61        else {
62            $trouble = 1;
63            die qq(Could not handle $ARGV[0]);
64      }      }
65  }  }
66    
67    my $taxon_ID;
68    my $org_dir;
69  (  (
70   ($genome = shift(@ARGV)) &&   ($taxon_ID = shift(@ARGV)) &&
71   ($dir    = shift(@ARGV))   ($org_dir    = shift(@ARGV))
72  )  )
73      || die $usage;      || die $usage;
74    
75  $prefixP = "fig|$genome.peg.";  $prefixP = qq(fig|$taxon_ID.peg.);
76  $prefixR = "fig|$genome.rna.";  $prefixR = qq(fig|$taxon_ID.rna.);
77    
78    #...Create paths down to PEG and RNA dirs
79    &FIG::verify_dir("$org_dir/Features/peg");
80    &FIG::verify_dir("$org_dir/Features/rna");
81    
82    my $fh_contigs;
83    my $contigs_file = qq($org_dir/contigs);
84    open($fh_contigs,  qq(>$contigs_file))
85        || die qq(could not open $contigs_file);
86    
87    my $fh_assigned_funcs;
88    my $assigned_funcs_file = qq($org_dir/assigned_functions);
89    open($fh_assigned_funcs,  qq(>$assigned_funcs_file))
90        || die qq(could not write-open $assigned_funcs_file);
91    
92    my $fh_annotations;
93    my $annotations_file = qq($org_dir/annotations);
94    open($fh_annotations,  qq(>$annotations_file))
95        || die qq(could not write-open $annotations_file);
96    
97    my $fh_ec_nums;
98    my $ec_nums_file = qq($org_dir/EC_numbers);
99    open($fh_ec_nums, qq(>$ec_nums_file))
100        || die" Could not write-open $ec_nums_file";
101    
102    my $fh_peg_tbl;
103    my $peg_tbl_file = qq($org_dir/Features/peg/tbl);
104    open($fh_peg_tbl,  qq(>$peg_tbl_file))
105        || die qq(could not open $peg_tbl_file);
106    
107    my $fh_peg_fasta;
108    my $peg_fasta_file = qq($org_dir/Features/peg/fasta);
109    open($fh_peg_fasta, qq(>$peg_fasta_file))
110        || die qq(could not open $peg_fasta_file);
111    
112    my $fh_rna_tbl;
113    my $rna_tbl_file = qq($org_dir/Features/rna/tbl);
114    open($fh_rna_tbl,      qq(>$rna_tbl_file))
115        || die qq(could not write-open $rna_tbl_file);
116    
117    my $fh_rna_fasta;
118    my $rna_fasta_file = qq($org_dir/Features/rna/fasta);
119    open($fh_rna_fasta,    qq(>$rna_fasta_file))
120        || die qq(could not write-open $rna_fasta_file);
121    
122    
123  &FIG::verify_dir("$dir/Features/peg");  open( TAXONOMY,  qq(>$org_dir/TAXONOMY))  || die qq(could not open $org_dir/TAXONOMY);
124  &FIG::verify_dir("$dir/Features/rna");  open( GENOME,    qq(>$org_dir/GENOME))    || die qq(could not open $org_dir/GENOME);
125  open(CONTIGS,">$dir/contigs") || die "could not open $dir/contigs";  
126  open(TBLPEG,">$dir/Features/peg/tbl") || die "could not open $dir/Features/peg/tbl";  open( PROJECT,   qq(>$org_dir/PROJECT))   || die qq(could not open $org_dir/PROJECT);
127  open(FASTAPEG,">$dir/Features/peg/fasta") || die "could not open $dir/Features/peg/fasta";  print PROJECT    qq($data_source\n);
 open(ASSIGNMENTS,">$dir/assigned_functions") || die "could not open $dir/assigned_functions";  
 open(TBLRNA,">$dir/Features/rna/tbl") || die "could not open $dir/Features/rna/tbl";  
 open(FASTARNA,">$dir/Features/rna/fasta") || die "could not open $dir/Features/rna/fasta";  
   
 open(TAX,">$dir/TAXONOMY") || die "could not open $dir/TAXONOMY";  
 open(GENOME,">$dir/GENOME") || die "could not open $dir/GENOME";  
 open(PROJECT,">$dir/PROJECT") || die "could not open $dir/PROJECT";  
 print PROJECT "Parsed NCBI Reference Sequences\n";  
128  close(PROJECT);  close(PROJECT);
129    
130  $idNp = 1;  my $idNp = 1;
131  $idNr = 1;  my $idNr = 1;
132    
133  $/ = "\n//\n";  $/ = qq(\n//\n);
134    my $record;
135    my $contigs = {};
136  while (defined($record = <STDIN>)) {  while (defined($record = <STDIN>)) {
137      $record =~ s/^\s+//os;      $record =~ s/^\s+//os;
138      last unless $record;      last unless $record;
139    
140      if ($record !~ m/LOCUS\s+(\S+)/os) {      if ($record !~ m/LOCUS\s+(\S+)/os) {
141          print STDERR "No LOCUS line for record begining with:\n"          print STDERR qq(No LOCUS line for record begining with:\n)
142              , substr($record, 0, 160), "\n\n";              , substr($record, 0, 160), qq(\n\n);
143      }      }
144      else {      else {
145          $id = $1;          my $id = $1;
146          if ($record =~ m/\nORIGIN\b(.*?)(\/\/|LOCUS)/os) {          if ($record =~ m/\nORIGIN\b(.*?)(\/\/|LOCUS)/os) {
147              $seq = $1;              my $seq = $1;
148              $seq =~ s/\s//ogs;              $seq =~ s/\s//ogs;
149              $seq =~ s/\d//og;              $seq =~ s/\d//og;
150              undef $contigs;              undef $contigs;
151              $contigs->{$id} = $seq;              $contigs->{$id} = $seq;
152              &FIG::display_id_and_seq($id,\$seq,\*CONTIGS);              &FIG::display_id_and_seq($id, \$seq, $fh_contigs);
153          }          }
154          else {          else {
155              warn "could not find contig sequence for $id\n";              warn qq(could not find contig sequence for $id\n);
156              next;              next;
157          }          }
158    
159            my ($taxon_ID, $taxonomy, $written_genome);
160          if ($record =~ /\n {0,4}ORGANISM\s+(\S[^\n]+(\n\s{10,14}\S[^\n]+)*)/os) {          if ($record =~ /\n {0,4}ORGANISM\s+(\S[^\n]+(\n\s{10,14}\S[^\n]+)*)/os) {
161              my $block = $1;              my $block = $1;
162              my @lines = split(/\n/,$block);              my @lines = split(/\n/,$block);
# Line 113  Line 172 
172                  $i++;                  $i++;
173              }              }
174    
175              $genome = join(" ",map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @genome);              $genome = join(qq( ),map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @genome);
176              $tax    = join(" ",map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @full_tax);              $taxonomy    = join(qq( ),map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @full_tax);
177    
178              $tax =~ s/\n\s+//og;              $taxonomy =~ s/\n\s+//og;
179              $tax =~ s/ {2,}/ /og;              $taxonomy =~ s/ {2,}/ /og;
180              $tax =~ s/\.$//o;              $taxonomy =~ s/\.$//o;
181              $tax = $tax . "; $genome";              $taxonomy = $taxonomy . qq(; $genome);
182    
183              if (! $written_genome) {              if (! $written_genome) {
184                  if ($force_genome) { $genome = $force_genome; }                  if ($force_bioname) { $taxon_ID = $force_bioname; }
185                  print GENOME "$genome\n";                  print GENOME qq($taxon_ID\n);
186                  close(GENOME);                  close(GENOME);
187    
188                  if ($force_tax) { $tax = $force_tax; }                  if ($force_taxonomy) { $taxonomy = $force_taxonomy; }
189                  print TAX "$tax\n";                  print TAXONOMY qq($taxonomy\n);
190                  close(TAX);                  close(TAXONOMY);
191    
192                  $written_genome = "$tax\t$genome";                  $written_genome = qq($taxonomy\t$bioname);
193              }              }
194              elsif (($written_genome ne "$tax\t$genome") && (!$force_genome || !$force_tax)) {              elsif (($written_genome ne qq($taxonomy\t$bioname)) && (!$force_bioname || !$force_taxonomy)) {
195                  print STDERR "serious mismatch in GENOME/TAX for $id\n$written_genome\n$tax\t$genome\n\n";                  print STDERR qq(serious mismatch in GENOME/TAX for $id\n$written_genome\n$taxonomy\t$bioname\n\n);
196              }              }
197          }          }
198    
199          while ($record =~ m/\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/ogs) {          while ($record =~ m/\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/ogs) {
200              $cds = $1;              my $cds = $1;
201              if (($cds !~ m/\/pseudo/o) &&              if (($cds !~ m/\/pseudo/o) &&
202                  (($cds !~ m/\/exception/o) || ($cds =~ m/\/translation/o))                  (($cds !~ m/\/exception/o) || ($cds =~ m/\/translation/o))
203                  ) {                  ) {
204                  &process_cds($id,\$cds,$prefixP,\$idNp,$contigs,\*TBLPEG,\*FASTAPEG,\*ASSIGNENTS);                  &process_cds($id, \$cds, $prefixP, \$idNp, $contigs, $fh_peg_tbl, $fh_peg_fasta, $fh_assigned_funcs, $fh_annotations, $fh_ec_nums);
205              }              }
206          }          }
207    
208          while ($record =~ m/\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/ogs)          while ($record =~ m/\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/ogs)
209          {          {
210              $rna = $1;              $rna = $1;
211              &process_rna($id,\$rna,$prefixR,\$idNr,$contigs,\*TBLRNA,\*FASTARNA);              &process_rna($id, \$rna, $prefixR, \$idNr, $contigs, $fh_rna_tbl, $fh_rna_fasta);
212            }
213        }
214    }
215    close($fh_contigs);
216    close($fh_peg_tbl);
217    close($fh_peg_fasta);
218    close($fh_rna_tbl);
219    close($fh_rna_fasta);
220    close($fh_assigned_funcs);
221    close($fh_annotations);
222    close($fh_ec_nums);
223    
224    if (!-s $assigned_funcs_file)  {
225        warn qq(no assigned_functions in $org_dir\n);
226        unlink("$assigned_funcs_file") unless $ENV{DEBUG};
227    }
228    
229    if (!-s $ec_nums_file)  {
230        warn qq(no EC_numbers in $org_dir\n);
231        unlink("$ec_nums_file") unless $ENV{DEBUG};
232    }
233    
234    if ((!-s $peg_tbl_file) || (!-s $peg_fasta_file)) {
235        $trouble = 1;
236        warn qq(ERROR: no PEGs in $org_dir\n);
237        system(qq(rm -rf $org_dir/Features/peg)) unless $ENV{DEBUG};
238    }
239    
240    if ((!-s $rna_tbl_file) || (!-s $rna_fasta_file)) {
241        warn qq(WARNING: no RNAs in $org_dir\n);
242        system(qq(rm -rf $org_dir/Features/rna)) unless $ENV{DEBUG};;
243          }          }
244    if (!-s $contigs_file) {
245        $trouble = 1;
246        unlink($contigs_file);
247        print STDERR qq(no contigs in $org_dir\n);
248      }      }
249    
250    if ($trouble) {
251        warn qq(Genome directory $org_dir is corrupt --- deleting\n\n);
252        system(qq(rm -rf $org_dir)) unless $ENV{DEBUG};
253  }  }
 close(CONTIGS);  
 close(TBLPEG);  
 close(FASTAPEG);  
 close(ASSIGNMENTS);  
 close(TBLRNA);  
 close(FASTARNA);  
   
 if (! -s "$dir/contigs")            { unlink("$dir/contigs"); print STDERR "no contigs in $dir\n"; }  
 if (! -s "$dir/assigned_functions") { unlink("$dir/assigned_functions"); print STDERR "no assigned_functions in $dir\n"; }  
 if (! -s "$dir/Features/peg/tbl")   { system "rm -rf $dir/Features/peg"; print STDERR "no PEGs in $dir\n"; }  
 if (! -s "$dir/Features/rna/tbl")   { system "rm -rf $dir/Features/rna"; print STDERR "no RNAs in $dir\n"; }  
 if ((! -s "$dir/contigs") && (! -s "$dir/Features/peg/tbl")) { system "rm -rf $dir" }  
254    
255    exit($trouble);
256    
257  sub process_cds {  sub process_cds {
258      my($contigID,$cdsP,$prefix,$idNp,$contigs,$fh_tbl,$fh_fasta,$fh_ass) = @_;      my ($contigID, $cdsP, $prefix, $idNp, $contigs, $fh_tbl, $fh_fasta, $fh_ass, $fh_annot, $fh_ec_nums) = @_;
259      my($loc,@aliases,$func,$trans,$id,$precise,$dna,$prot);      my ($id, $prot);
260    
261      ++$recnum;      ++$recnum;
262      ($loc,$precise)  = &get_loc($contigID,$cdsP);      my ($loc, $precise)  = &get_loc($contigID,$cdsP);
263      @aliases = &get_aliases($cdsP);      my @aliases = &get_aliases($cdsP);
264      $func    = &get_func($cdsP);      my @ec_nums = &get_ec_numbers($cdsP);
265      $trans   = &get_trans($cdsP);      my ($func, $notes) = &get_func($cdsP);
266        my $trans   = &get_trans($cdsP);
267    
268      if ($loc || $trans) {      if ($loc || $trans) {
269          if ($dna  = &FIG::extract_seq($contigs,$loc)) {          my $dna = &FIG::extract_seq($contigs,$loc);
270            if ($dna) {
271              if ($precise) {              if ($precise) {
272                  $prot = &FIG::translate($dna,undef,1);                  if (not $trans) {
273                  if ($trans && (uc $prot ne uc $trans)) {                      warn qq(WARNING: Translation missing for CDS $recnum; attempting to generate translation\n);
274  #               print STDERR "BAD TRANSLATION: $contigID $$cdsP\n";  
275  #               print STDERR &Dumper($trans,$prot,$dna),"\n";                      my $genetic_code = 11;
276                        if ($$cdsP =~ m/\/transl_table=(\d+)/o) {
277                            $genetic_code = $1;
278                  }                  }
279    
280                        if (not defined($trans_table{$genetic_code})) {
281                            $trans_table{$genetic_code} = &FIG::genetic_code($genetic_code);
282              }              }
283    
284                        my $prot = &FIG::translate($dna, $trans_table{$genetic_code}, 1);
285    
286                        if ($prot !~ m/\*/o) {
287                            $trans = $prot;
288          }          }
289          else {          else {
290              warn "Could not get DNA sequence for CDS at $loc for entry:\n$$cdsP\nof record begining with:\n"                          warn qq(Translation contains STOPs, changing to \'x\'s, for CDS $recnum:\n$$cdsP\n);
291                  , substr($record, 0, 160), "\n\n";                      }
292                    }
293                }
294            }
295            else {
296                warn qq(Could not get DNA sequence for CDS at $loc for entry:\n$$cdsP\nof record begining with:\n)
297                    , substr($record, 0, 160), qq(\n\n);
298          }          }
299    
300          $id = $prefix . "$$idNp";          $id = $prefix . qq($$idNp);
301          $$idNp++;          ++$$idNp;
302          print $fh_tbl "$id\t$loc\t",join("\t",@aliases),"\n";          print $fh_tbl qq($id\t$loc\t",join("\t",@aliases),"\n);
303    
304          if ($func) {          if ($func) {
305              print $fh_ass "$id\t$func\n";              print $fh_ass qq($id\t$func\n);
306            }
307    
308            if (@ec_nums) {
309                print $fh_ec_nums (join(qq(\t), ($id, @ec_nums)), qq(\n));
310            }
311    
312            foreach my $note (@$notes) {
313                &make_annotation($id, $note, $fh_annot);
314          }          }
315    
316          if ($trans) {          if ($trans) {
317              &FIG::display_id_and_seq($id,\$trans,$fh_fasta);              &FIG::display_id_and_seq($id,\$trans,$fh_fasta);
318          }          }
319          else {          else {
320              warn "No translation for $id";              die qq(No translation for $id);
321          }          }
322      }      }
323      else {      else {
324          print &Dumper($cdsP); die "aborted";          print &Dumper($cdsP); die qq(aborted);
325      }      }
326  }  }
327    
# Line 256  Line 369 
369          }          }
370          else {          else {
371              print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);              print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
372              die "aborted";              die qq(aborted);
373          }          }
374    
375          my @locs = split(m/,/o, $loc);          my @locs = split(m/,/o, $loc);
# Line 265  Line 378 
378          $loc = join(",",@locs);          $loc = join(",",@locs);
379      }      }
380      else {      else {
381          print STDERR &Dumper($cdsP); die "could not parse location";          print STDERR &Dumper($cdsP); die qq(could not parse location);
382          die "aborted";          die qq(aborted);
383      }      }
384      return ($loc,$precise);      return ($loc,$precise);
385  }  }
# Line 293  Line 406 
406              }              }
407          }          }
408          else {          else {
409              die "could not parse $locs->[$n]";              die qq(could not parse $locs->[$n]);
410          }          }
411      }      }
412  }  }
# Line 302  Line 415 
415  sub conv {  sub conv {
416      my($pieces,$n,$contigID) = @_;      my($pieces,$n,$contigID) = @_;
417    
418      if ($pieces->[$n]->[0] eq "loc") {      if ($pieces->[$n]->[0] eq qq(loc)) {
419          return join("_",$contigID,@{$pieces->[$n]}[1..2]);          return join("_",$contigID,@{$pieces->[$n]}[1..2]);
420      }      }
421      elsif ($pieces->[$n]->[0] eq "join") {      elsif ($pieces->[$n]->[0] eq qq(join)) {
422          return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);          return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
423      }      }
424      elsif ($pieces->[$n]->[0] eq "complement") {      elsif ($pieces->[$n]->[0] eq qq(complement)) {
425          return join(",",&complement(join(",", map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;          return join(",",&complement(join(",", map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
426      }      }
427  }  }
# Line 323  Line 436 
436              $loc = join("_",($1,$3,$2));              $loc = join("_",($1,$3,$2));
437          }          }
438          else {          else {
439              confess "Bad location: $loc";              confess qq(Bad location: $loc);
440          }          }
441      }      }
442      return join(",",@locs);      return join(",",@locs);
# Line 331  Line 444 
444    
445  sub get_aliases {  sub get_aliases {
446      my($cdsP) = @_;      my($cdsP) = @_;
447      my($db_ref);      my($id, $prefix, $alias, $db_ref);
448    
449      my @aliases = ();      my @aliases = ();
450      while ($$cdsP =~ m/\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/ogs) {      while ($$cdsP =~ m/\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/ogs) {
451        my $id;          ($type, $alias) = ($1, $alias);
452    
453        # define prefixes for different types of ids        # define prefixes for different types of ids
454            if ($type eq qq(locus_tag)){
455        if ($1 eq "locus_tag"){              $id = qq(locus|$alias);
         $id = "locus|$2";  
456        }        }
457        elsif ( $1 eq "protein_id" ) {          elsif ( $type eq qq(protein_id) ) {
458          $id = "protein_id|$2";              $id = qq(protein_id|$alias);
459        }        }
460        elsif  ( $1 eq "gene" ){          elsif  ( $type eq qq(gene) ){
461          $id = "gene_name|$2";              $id = qq(gene_name|$alias);
462        }        }
463        else{        else{
464          $id = $2;              $id = $alias;
465        }        }
466    
   
467        push(@aliases,$id);        push(@aliases,$id);
468      }      }
469    
# Line 367  Line 478 
478      return @aliases;      return @aliases;
479  }  }
480    
481    sub get_ec_numbers {
482        my ($cdsP) = @_;
483        my @ec_numbers = ($$cdsP =~ m{/EC_number=\"([0-9]+\.[0-9\-]+\.[0-9\-]+\.[0-9\-]+)\"}ogs);
484        return @ec_numbers;
485    }
486    
487  sub get_trans {  sub get_trans {
488      my($cdsP) = @_;      my($cdsP) = @_;
489      my $tran = "";      my $tran = qq();
490    
491      if ($$cdsP =~ m/\/translation=\"([^"]*)\"/os) {      if ($$cdsP =~ m/\/translation=\"([^"]*)\"/os) {
492          $tran = $1;          $tran = $1;
# Line 384  Line 501 
501    
502  sub get_func {  sub get_func {
503      my($cdsP) = @_;      my($cdsP) = @_;
     my $func = "";  
504    
505      print STDERR "\nRecord $recnum:\n$$cdsP\n" if $ENV{VERBOSE};      my $functions  = [];
506      if (($$cdsP =~ m/\/function=\"([^"]*)\"/os) && ($func = $1) && &ok_func($func)) {      my $products   = [];
507          print STDERR "Branch 1: $func\n" if $ENV{VERBOSE};      my $prot_descs = [];
508      }      my $notes      = [];
509      elsif (($$cdsP =~ m/\/product=\"([^"]*)\"/os)  
510          && ($func = $1) && &ok_func($func)      print STDERR qq(\nRecord $recnum:\n$$cdsP\n) if $ENV{VERBOSE};
511          && (($func =~ m/ /o) || ($$cdsP !~ m/\/note/o))  
512          ) {      @$functions = ($$cdsP =~ m/\/function=\"([^\"]*)\"/ogs);
513          print STDERR "Branch 2: $func\n" if $ENV{VERBOSE};      if ($ENV{VERBOSE} && @$functions) {
514            print STDERR (qq(Functions: ), ((@$functions > 1) ? qq(\n) : qq()));
515            print STDERR (join(qq(\n), @$functions), qq(\n));
516      }      }
517      elsif (($$cdsP =~ m/\/note=\"([^"]*)\"/os) && ($func = $1) && &ok_func($func)) {  
518          print STDERR "Branch 3: $func\n" if $ENV{VERBOSE};      @$products = ($$cdsP =~ m/\/product=\"([^\"]*)\"/ogs);
519        if ($ENV{VERBOSE} && @$products) {
520            print STDERR (qq(Products: ), ((@$products > 1) ? qq(\n) : qq()));
521            print STDERR (join(qq(\n), @$products), qq(\n));
522      }      }
523      else {  
524          print STDERR "No non-hypo found\n" if $ENV{VERBOSE};      @$prot_descs = ($$cdsP =~ m/\/prot_desc=\"([^\"]*)\"/ogs);
525        if ($ENV{VERBOSE} && @$prot_descs) {
526            print STDERR (qq(Prot_Descs: ), ((@$prot_descs > 1) ? qq(\n) : qq()));
527            print STDERR (join(qq(\n), @$prot_descs), qq(\n));
528        }
529    
530        @$notes = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/note=\"([^\"]*)\"/ogs);
531        if ($ENV{VERBOSE} && @$notes) {
532            print STDERR (qq(Notes: ), ((@$notes > 1) ? qq(\n) : qq()));
533            print STDERR (join(qq(\n), @$notes), qq(\n));
534      }      }
     $func =~ s/\s+/ /ogs;  
     print STDERR     "      --> $func\n" if $ENV{VERBOSE};  
535    
536      $func = &fixup_func($func);      @$products = grep { !m/^hypothetical\s+protein$/io } @$products;
     print STDERR     "Retval:   $func\n\n" if $ENV{VERBOSE};  
537    
538      return $func;      my $func  = qq();
539        my $annotations = [];
540        if (@$products) {
541            $func = join(q( / ), @$products);
542            @$annotations = (@$functions, @$prot_descs, @$notes);
543        }
544        elsif (@$functions) {
545            $func =join(q( / ), @$functions);
546            @$annotations =(@$prot_descs, @$notes);
547        }
548        elsif (@$prot_descs) {
549            $func =join(q( / ), @$prot_descs);
550            @$annotations = @$notes;
551        }
552        else {
553            $func = qq(hypothetical protein);
554            @$annotations = @$notes;
555        }
556    
557        return ($func, $annotations);
558  }  }
559    
560  sub fixup_func {  sub fixup_func {
# Line 452  Line 598 
598    
599      if ($loc) {      if ($loc) {
600          if ($dna  = &FIG::extract_seq($contigs,$loc)) {          if ($dna  = &FIG::extract_seq($contigs,$loc)) {
601              $id = $prefix . "$$idNr";              $id = $prefix . qq($$idNr);
602              $$idNr++;              $$idNr++;
603    
604              if (! $func ) {              if (! $func ) {
605                  warn "WARNING: could not get func $$cdsP\n" if $ENV{VERBOSE};                  warn qq(WARNING: could not get func $$cdsP\n) if $ENV{VERBOSE};
606              }              }
607              else {              else {
608                  print $fh_tbl "$id\t$loc\t$func\n";                  print $fh_tbl qq($id\t$loc\t$func\n);
609                  &FIG::display_id_and_seq($id,\$dna,$fh_dna);                  &FIG::display_id_and_seq($id,\$dna,$fh_dna);
610              }              }
611          }          }
612          else {          else {
613              warn "Could not get DNA sequence for RNA at $loc for entry\n$$cdsP\nof record begining with:\n"              warn qq(Could not get DNA sequence for RNA at $loc for entry\n$$cdsP\nof record begining with:\n)
614              , substr($record, 0, 160), "\n\n";              , substr($record, 0, 160), qq(\n\n);
615          }          }
616      }      }
617      else {      else {
618          print &Dumper($cdsP); die "aborted";          print &Dumper($cdsP); die qq(aborted);
619      }      }
620  }  }
621    
622  sub get_rna_func {  sub get_rna_func {
623      my($cdsP) = @_;      my($cdsP) = @_;
624      my $func = "";      my $func = qq();
625    
626      if ($$cdsP =~ m/\/product=\"([^"]*)\"/os) {      if ($$cdsP =~ m/\/product=\"([^"]*)\"/os) {
627          $func = $1;          $func = $1;
# Line 491  Line 637 
637      }      }
638      return $func;      return $func;
639  }  }
640    
641    sub make_annotation {
642        my ($fid, $note, $fh_ann) = @_;
643    #   print STDERR Dumper($fid, $note, $fh_ann);
644    
645        print $fh_ann ($fid, qq(\n));
646        print $fh_ann (time, qq(\n));
647        print $fh_ann qq(parse_genbank\n);
648        print $fh_ann ($note, qq(\n));
649        print $fh_ann qq(//\n);
650    
651        return 1;
652    }

Legend:
Removed from v.1.17  
changed lines
  Added in v.1.18

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3