#!/usr/bin/perl -w

# By HongLu @ 2011/02/05
# Email: luhong@iastate.edu

# Retrieve protein sequences from GenBank files

use strict;
use Getopt::Std;

my %opt = ();
getopts("g:o:", \%opt);
my $usage = "perl $0
    <-g> GenBank file
    [-o] outfile, the default name is {GenBank file}.faa\n";

my ($GenBank, $outfile);
if (defined $opt{g}) { $GenBank = $opt{g}; } else { print $usage; exit; }
if (defined $opt{o}) { $outfile = $opt{o}; } else { $outfile = $GenBank.".faa"; }

my ($step, $guard, $space) = (0, 0, "");
my ($geneName, $geneCoord, $geneID, $proSeq, $function) = ("", "", "", "", "");
# $step
# 0: initial; 1: gene; 2: CDS; 3: end
open(IN, $GenBank)||die ":$!";
open(OUT,">$outfile")||die ":$!";
while (<IN>) {
  chomp;
  if (/^(\s+)source/) {
	if ($space eq "") {
	  $space = $1;
	}
  }
  if (/^ORIGIN/) {
    unless ($geneName eq "") { # Not the first one
      $geneID = ($geneID eq "")?"xxxN":$geneID;
      if ($step eq "2") {
        printSeq($geneName, $geneCoord, $geneID, $proSeq, $function);
      }
      ($geneName, $geneCoord, $geneID, $proSeq, $function) = ("", "", "", "", "");
    }
    $step = 3;
    last;
  }
  if (/^${space}gene\s+(\S+)/) {
    unless ($geneName eq "") { # Not the first one
      $geneID = ($geneID eq "")?"xxxN":$geneID;
      if ($step eq "2") {
        printSeq($geneName, $geneCoord, $geneID, $proSeq, $function);
      }
      ($geneName, $geneCoord, $geneID, $proSeq, $function) = ("", "", "", "", "");
    }
    ($step, $geneCoord) = (1, $1);
    next;
  }
  if (($step eq "1")&&((/^\s+\/locus_tag=\"(\S+)\"/)||(/^\s+\/locus_tag=\"(\S+)\;\s*\S+\"/))) {
    if ($geneName ne "") {
      next;
    } else {
      $geneName = $1;
    }
    next;
  }
  if (($step eq "1")&&(/^${space}CDS\s+\S+/)) {
    $step = 2;
    next;
  }
  if ((($step eq 2)||($step eq 1))&&(/^\s+\/product\=\"(.*)/)) {
    if ($function ne "") {
      next;
    }
    $function = $1;
    $guard    = 0;
    if ($function=~/\"$/) {
      chop($function);
    } else {
      $guard = 1;
    }
    next;
  }
  if ($guard eq "1") {
    if (/\"$/) {
      $function = $function.$`;
      $guard    = 0;
      $function =~s/\s+/ /g;
    } else {
      $function = $function.$_;
    }
    next;
  }
  if (($step eq "2")&&(/^\s+\/gene\=\"(.*)\"/)) {
    $geneID = $1;
    next;
  }
  if (($step eq 2)&&(/^\s+\/translation\=\"(\S+)/)) {
    $proSeq = $1;
    $guard  = 0;
    if ($proSeq=~/\"$/) {
      chop($proSeq);
    } else {
      $guard = 2;
    }
    next;
  }
  if ($guard eq "2") {
    if (/\"$/) {
      $proSeq = $proSeq.$`;
      $proSeq =~s/\s+//g;
      $proSeq =~s/\"//g;
      $guard  = 0;
    } else {
      $proSeq = $proSeq.$_;
    }
    next;
  }
}
close(IN)||die ":$!";
close(OUT)||die ":$!";

#***********************************************************
sub printSeq {
  my ($geneName, $geneCoord, $geneID, $proSeq, $function) = @_;
  print OUT ">$geneName    $geneID    $function\n";
  print OUT $proSeq."\n";
}
