prokka-genbank_to_fasta_db

#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
use Data::Dumper;

my(@Options, $verbose, $format, $gcode, $hypo, $idtag, $sep, $blank, $pseudo, $minlen);
setOptions();
($gcode < 1 or $gcode > 24) and die "Invalid genetic code, must be 1..24\n";
print STDERR "Using genetic code table $gcode.\n";

my $in = Bio::SeqIO->new(-fh=>\*ARGV, -format=>$format);
my $out = Bio::SeqIO->new(-fh=>\*STDOUT, -format=>'Fasta');

while (my $seq = $in->next_seq) {
  print STDERR "\rParsing: ",$seq->display_id;
  my $counter = 0;
  for my $f ($seq->get_SeqFeatures)   {

    next unless $f->primary_tag eq 'CDS';

    next unless $f->length >= $minlen;

    next if !$pseudo and $f->has_tag('pseudo');

    my $prod = TAG($f, 'product') || $blank;

    next if !$hypo and $prod eq 'hypothetical protein';

    my $cds = $f->spliced_seq;  # don't forget eukaryotes!

    my $id = TAG($f, $idtag) or die "feature #$counter does not have /$idtag tag!";
    $counter++;
    print STDERR "Processing: [$counter] ", $seq->display_id, " | $id | $prod\n" if $verbose;

    # HANDLE CODON START FOR FUZZY FEATURES!
    if ($f->has_tag('codon_start')) {
      my($cs) = $f->get_tag_values('codon_start');
      if ($cs != 1) {
        print STDERR "/codon_start=$cs - trimming mRNA!";
        $cds = $cds->trunc($cs, $cds->length);
      }
    }
    #END


    # DNA -> AA
    $cds = $cds->translate(-codontable_id=>$gcode, -complete => 1);

    my $ec = TAG($f, 'EC_number') || $blank;
    my $gene = TAG($f, 'gene') || $blank;

    $cds->desc("$ec$sep$gene$sep$prod");
    $cds->display_id($id);

    $out->write_seq($cds);
  }
  print STDERR "\n";
}
print STDERR "\nDone\n";

#----------------------------------------------------------------------

sub TAG {
  my($f, $tag) = @_;
  return unless $f->has_tag($tag);
  # Seems new submissions have 2 protein IDs - one useful and one annoying!
  # /protein_id="REF_PRJNA193299:EFAU085_p1001"
  # /protein_id="YP_008390691.1"
  return join(';', $f->get_tag_values($tag)) if ($tag eq 'EC_number'); # concatenate EC_numbers via ';' if several present
  return ( grep { ! m/PRJNA/ } $f->get_tag_values($tag) ) [0];
}

#----------------------------------------------------------------------
# Option setting routines

sub setOptions {
  use Getopt::Long;

  @Options = (
    {OPT=>"help",      VAR=>\&usage,               DESC=>"This help"},
    {OPT=>"verbose!",  VAR=>\$verbose, DEFAULT=>0, DESC=>"Verbose progress"},
    {OPT=>"format=s",  VAR=>\$format,  DEFAULT=>'genbank', DESC=>"Input format"},
    {OPT=>"gcode=i",  VAR=>\$gcode, DEFAULT=>1, DESC=>"Genetic code / Translation table (1..24)"},
    {OPT=>"idtag=s",   VAR=>\$idtag,   DEFAULT=>'protein_id', DESC=>"What tag to use as Fasta ID"},
    {OPT=>"sep=s",     VAR=>\$sep,     DEFAULT=>'~~~', DESC=>"Separator between EC/gene/product" },
    {OPT=>"blank=s",   VAR=>\$blank,   DEFAULT=>'',    DESC=>"Replace empty EC/gene/product with this"},
    {OPT=>"pseudo!",   VAR=>\$pseudo,  DEFAULT=>0,     DESC=>"Include /pseudo genes"},
    {OPT=>"hypo!",     VAR=>\$hypo,    DEFAULT=>0,     DESC=>"Include 'hypothetical protein' genes"},
    {OPT=>"minlen=i",  VAR=>\$minlen,  DEFAULT=>0,     DESC=>"Minimum peptide length"},
  );

  (!@ARGV) && (usage());

  &GetOptions(map {$_->{OPT}, $_->{VAR}} @Options) || usage();

  # Now setup default values.
  foreach (@Options) {
    if (defined($_->{DEFAULT}) && !defined(${$_->{VAR}})) {
      ${$_->{VAR}} = $_->{DEFAULT};
    }
  }
}

sub usage {
  print "Usage: $0 [options] [genome1.gbk ...] > proteins.faa\n";
  foreach (@Options) {
    printf "  --%-13s %s%s.\n",$_->{OPT},$_->{DESC},
           defined($_->{DEFAULT}) ? " (default '$_->{DEFAULT}')" : "";
  }
  exit(1);
}

#----------------------------------------------------------------------