#!/usr/bin/perl
use strict;
use warnings;
##
## seq_analysis.pl
##
##   performing several analyses on input DNA sequences
##   and providing all-frames translation
##

## Variables Declaration and Initialization

our ($dna_id, $dna_seq, $dna_revseq, $dna_seqlen);


## Loading Sequence

while (<>) {

	next if /^\#/o;
	next if /^\s*$/o;

	chomp;

	my $line = $_;

	$line =~ /^\>([^\s]+)\s*/o && do {

		$dna_id = defined($1) ? $1 : 'unknown_dna_seq';
		$dna_seq = '';
		next;

	};

	$dna_seq .= $line;

}; # while STDIN


## Basic Sequence Features

our %NUCLEOTIDES = ();

$dna_seqlen = length($dna_seq);

for (my $i = 0; $i < $dna_seqlen; $i++) { # Get Nucleotide Counts
    my $char = uc(substr($dna_seq,$i,1));
    $NUCLEOTIDES{$char}++;
};

my $gc_content = (($NUCLEOTIDES{G} + $NUCLEOTIDES{C}) / $dna_seqlen) * 100;

print  STDOUT '=' x 80,"\n\n";
printf STDOUT "SEQUENCE: %s\n\n".
              "--> STATS:\n".
              "  Length: %dbp\n".
              "     G+C: %5.2f%%\n\n",
	          $dna_id, $dna_seqlen, $gc_content;


## Reverse Complement

( $dna_revseq = reverse $dna_seq ) =~ tr/acgtACGT/tgcaTGCA/;


## Translate Sequence

our %CODON_TABLE = ();

while (<DATA>) { # Loading Translation Hash

    next if $_ =~ m/^\#/o;

	my @F = split /\s+/o, $_;

	for (my $i = 0; $i < scalar(@F); $i+=2) {

		$CODON_TABLE{$F[$i]} = $F[($i + 1)];

	}; # for $i

}; # while DATA


foreach my $fr (0, 1) { # Do forward and reverse

	my $seqref = ($fr == 0) ? \$dna_seq : \$dna_revseq;

	printf STDOUT "--> %s Translations:\n",
                  (($fr == 0) ? "FORWARD" : "REVERSE");

	for (my $phase = 0; $phase < 3; $phase++) {

		my $protein_seq = '';

		for (my $i = $phase; $i < $dna_seqlen; $i += 3) {

			my $codon = substr($$seqref,$i,3);

			last if length($codon) < 3;

			if (defined($CODON_TABLE{$codon})) {

				$protein_seq .= $CODON_TABLE{$codon};

			} else {

				$protein_seq .= 'X'; # Any wildcard provided by
		                             # the corresponding annotation standard...

			}; # 

		}; # for $i

		printf STDOUT "   Phase %d : %s\n", $phase, $protein_seq;

	}; # for $phase

	print STDOUT "\n";

}; # foreach $refseq

print STDOUT '=' x 80,"\n\n";


## The End
exit(0); # optional, but to be polite with the parent shell...


## Anyting following the '__END__' or '__DATA__' tokens is ignored 

__END__
# Standard Codon Table -- shorter version
TTT F  TTC F  TTA L  TTG L  CTT L  CTC L  CTA L  CTG L
ATT I  ATC I  ATA I  ATG M  GTT V  GTC V  GTA V  GTG V
TCT S  TCC S  TCA S  TCG S  CCT P  CCC P  CCA P  CCG P
ACT T  ACC T  ACA T  ACG T  GCT A  GCC A  GCA A  GCG A
TAT Y  TAC Y  TAA *  TAG *  CAT H  CAC H  CAA Q  CAG Q
AAT N  AAC N  AAA K  AAG K  GAT D  GAC D  GAA E  GAG E
TGT C  TGC C  TGA *  TGG W  CGT R  CGC R  CGA R  CGG R
AGT S  AGC S  AGA R  AGG R  GGT G  GGC G  GGA G  GGG G
