[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.8, Sat Jun 12 11:15:56 2004 UTC revision 1.9, Fri Jul 16 21:54:14 2004 UTC
# Line 1  Line 1 
1  use FIG;  use FIG;
2    
3  $usage = "usage: parse_genbank Genome Dir < genbank.entry";  $usage = "usage: parse_genbank Genome Dir < genbank.entry";
4    while ($ARGV[0] =~ m/^-/)
5    {
6        if    ($ARGV[0] =~ m/^-{1,2}genome=\S+/)
7        {
8            $force_genome =  shift @ARGV;
9            $force_genome =~ s/^-{1,2}genome=\"?//;
10            $force_genome =~ s/\"$//;
11            $force_genome =~ s/\s+/ /g;
12    
13            print STDERR "\nGenome bioname will be taken as: \"$force_genome\"\n";
14        }
15        elsif ($ARGV[0] =~ m/^-{1,2}tax=\S+/)
16        {
17            $force_tax =  shift @ARGV;
18            $force_tax =~ s/^-{1,2}tax=\"?//;
19            $force_tax =~ s/\"$//;
20            $force_tax =~ s/\s+/ /g;
21            if ($force_tax !~ m/\.$/) { $force_tax .= "."; }
22    
23            print STDERR "\nTaxonomy will be taken as: \"$force_tax\"\n";
24        }
25        else
26        {
27            die "Could not handle $ARGV[0]";
28        }
29    }
30    
31  (  (
32   ($genome = shift(@ARGV)) &&   ($genome = shift(@ARGV)) &&
33   ($dir    = shift(@ARGV))   ($dir    = shift(@ARGV))
# Line 28  Line 55 
55  $idNp = 1;  $idNp = 1;
56  $idNr = 1;  $idNr = 1;
57    
58  $/ = "\nLOCUS";  $/ = "\n//\n";
59  while ($contig = <STDIN>)  while ($record = <STDIN>)
60  {  {
61      if ($contig =~ /\nACCESSION\s+(\S+)/s)      if ($record !~ /LOCUS\s+(\S+)/s)
62        {
63    
64            print STDERR "No LOCUS line for record begining with:\n"
65                , substr($record, 0, 160), "\n\n";
66        }
67        else
68      {      {
69          $id = $1;          $id = $1;
70          if ($contig =~ /ORIGIN(.*?)(\/\/|LOCUS)/s)          if ($record =~ /ORIGIN(.*?)(\/\/|LOCUS)/s)
71          {          {
72              $seq = $1;              $seq = $1;
73              $seq =~ s/\s//gs;              $seq =~ s/\s//gs;
# Line 45  Line 78 
78          }          }
79          else          else
80          {          {
81              warn "could not find sequence for $id\n";              warn "could not find contig sequence for $id\n";
82              next;              next;
83          }          }
84    
85          if ($contig =~ /\n {0,4}ORGANISM\s+(\S[^\n]+\S)((\n\s{10,14}\S[^\n]+\S)+)/s)          if ($record =~ /\n {0,4}ORGANISM\s+(\S[^\n]+\S)((\n\s{10,14}\S[^\n]+\S)+)/s)
86          {          {
87              $genome = $1;              $genome = $1;
88              $tax = $2;              $tax = $2;
# Line 57  Line 90 
90              $tax =~ s/ {2,}/ /g;              $tax =~ s/ {2,}/ /g;
91              if (! $written_genome)              if (! $written_genome)
92              {              {
93                    if ($force_genome) { $genome = $force_genome; }
94                  print GENOME "$genome\n";                  print GENOME "$genome\n";
95                  close(GENOME);                  close(GENOME);
96    
97                    if ($force_tax) { $tax = $force_tax; }
98                  print TAX "$tax\n";                  print TAX "$tax\n";
99                  close(TAX);                  close(TAX);
100    
101                  $written_genome = "$tax\t$genome";                  $written_genome = "$tax\t$genome";
102              }              }
103              elsif ($written_genome ne "$tax\t$genome")              elsif (($written_genome ne "$tax\t$genome") && (!$force_genome || !$force_tax))
104              {              {
105                  print STDERR "serious mismatch in GENOME/TAX for $id\n$written_genome\n$tax\t$genome\n\n";                  print STDERR "serious mismatch in GENOME/TAX for $id\n$written_genome\n$tax\t$genome\n\n";
106              }              }
107          }          }
108    
109          while ($contig =~ /\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/gs)          while ($record =~ /\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/gs)
110          {          {
111              $cds = $1;              $cds = $1;
112              if (($cds !~ /\/pseudo/) && (($cds !~ /\/exception/) || ($cds =~ /\/translation/)))              if (($cds !~ /\/pseudo/) && (($cds !~ /\/exception/) || ($cds =~ /\/translation/)))
# Line 79  Line 115 
115              }              }
116          }          }
117    
118          while ($contig =~ /\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/gs)          while ($record =~ /\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/gs)
119          {          {
120              $rna = $1;              $rna = $1;
121              &process_rna($id,\$rna,$prefixR,\$idNr,$contigs,\*TBLRNA,\*FASTARNA);              &process_rna($id,\$rna,$prefixR,\$idNr,$contigs,\*TBLRNA,\*FASTARNA);
# Line 102  Line 138 
138    
139  sub process_cds {  sub process_cds {
140      my($contigID,$cdsP,$prefix,$idNp,$contigs,$fh_tbl,$fh_fasta,$fh_ass) = @_;      my($contigID,$cdsP,$prefix,$idNp,$contigs,$fh_tbl,$fh_fasta,$fh_ass) = @_;
141      my($loc,@aliases,$func,$trans,$id,$precise);      my($loc,@aliases,$func,$trans,$id,$precise,$dna,$prot);
142    
143      ($loc,$precise)  = &get_loc($contigID,$cdsP);      ($loc,$precise)  = &get_loc($contigID,$cdsP);
144      @aliases = &get_aliases($cdsP);      @aliases = &get_aliases($cdsP);
# Line 111  Line 147 
147    
148      if ($loc || $trans)      if ($loc || $trans)
149      {      {
150          my $dna  = &FIG::extract_seq($contigs,$loc);          if ($dna  = &FIG::extract_seq($contigs,$loc))
151            {
152          if ($precise)          if ($precise)
153          {          {
154              my $prot = &FIG::translate($dna,undef,1);                  $prot = &FIG::translate($dna,undef,1);
155              if ($trans && (uc $prot ne uc $trans))              if ($trans && (uc $prot ne uc $trans))
156              {              {
157  #               print STDERR "BAD TRANSLATION: $contigID $$cdsP\n";  #               print STDERR "BAD TRANSLATION: $contigID $$cdsP\n";
158  #               print STDERR &Dumper($trans,$prot,$dna),"\n";  #               print STDERR &Dumper($trans,$prot,$dna),"\n";
159              }              }
160          }          }
161            }
162            else
163            {
164                warn "Could not get DNA sequence for CDS at $loc for entry:\n$$cdsP\nof record begining with:\n"
165                , substr($record, 0, 160), "\n\n";
166            }
167    
168          $id = $prefix . "$$idNp";          $id = $prefix . "$$idNp";
169          $$idNp++;          $$idNp++;
# Line 130  Line 173 
173          {          {
174              print $fh_ass "$id\t$func\n";              print $fh_ass "$id\t$func\n";
175          }          }
176    
177          if ($trans)          if ($trans)
178          {          {
179              &FIG::display_id_and_seq($id,\$trans,$fh_fasta);              &FIG::display_id_and_seq($id,\$trans,$fh_fasta);
180          }          }
181            else
182            {
183                warn "No translation for $id";
184            }
185      }      }
186      else      else
187      {      {
# Line 363  Line 411 
411    
412  sub process_rna {  sub process_rna {
413      my($contigID,$cdsP,$prefix,$idNr,$contigs,$fh_tbl,$fh_dna) = @_;      my($contigID,$cdsP,$prefix,$idNr,$contigs,$fh_tbl,$fh_dna) = @_;
414      my($loc,@aliases,$func,$trans,$id,$precise);      my($loc,@aliases,$func,$trans,$id,$precise,$DNA);
415    
416      ($loc,$precise)  = &get_loc($contigID,$cdsP);      ($loc,$precise)  = &get_loc($contigID,$cdsP);
417      $func    = &get_rna_func($cdsP);      $func    = &get_rna_func($cdsP);
418    
419      if ($loc)      if ($loc)
420      {      {
421          my $dna  = &FIG::extract_seq($contigs,$loc);          if ($dna  = &FIG::extract_seq($contigs,$loc))
422            {
423          $id = $prefix . "$$idNr";          $id = $prefix . "$$idNr";
424          $$idNr++;          $$idNr++;
425          if (! $func )          if (! $func )
# Line 385  Line 434 
434      }      }
435      else      else
436      {      {
437                warn "Could not get DNA sequence for RNA at $loc for entry\n$$cdsP\nof record begining with:\n"
438                , substr($record, 0, 160), "\n\n";
439            }
440        }
441        else
442        {
443          print &Dumper($cdsP); die "aborted";          print &Dumper($cdsP); die "aborted";
444      }      }
445  }  }

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3