#!/usr/bin/perl
use strict;
use warnings;
##
## seq_analysis_subs.pl
##
##   performing several analyses on input DNA sequences
##   and providing all-frames translation.
##
##   For this version, most of the code was moved into functions
##

## Variables Declaration and Initialization

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

&init_codon_table();


## Main Loop

($dna_id, $dna_seq) = &load_seq();

&spacer();

$dna_seqlen = &seq_stats($dna_id, \$dna_seq);

&revcomplement(\$dna_seq, \$dna_revseq);

&translate(\$dna_seq,    0); # 0 -> FORWARD | 1 -> REVERSE
&translate(\$dna_revseq, 1);

&spacer();

exit(0);


## Functions

sub spacer() { print STDOUT '=' x 80,"\n\n"; }

sub load_seq() { ## Loading Sequence

	my ($id, $seq);

	while (<>) {

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

		chomp;

		my $line = $_;

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

			$id = defined($1) ? $1 : 'unknown_dna_seq';
			$seq = '';
			next;

		};

		$seq .= $line;

	}; # while STDIN

	return ($id, $seq);

} # load_seq

sub seq_stats() { ## Basic Sequence Features

	my (%NUCLEOTIDES, $id, $rseq, $seqlen, $gc_content);

	($id, $rseq) = @_;

	%NUCLEOTIDES = ();

	$seqlen = length($$rseq);

	for (my $i = 0; $i < $seqlen; $i++) { # Get Nucleotide Counts

		my $char = uc(substr($$rseq, $i, 1));
		$NUCLEOTIDES{$char}++;

	};

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

	printf STDOUT "SEQUENCE: %s\n\n".
		          "--> STATS:\n".
		          "  Length: %dbp\n".
				  "     G+C: %5.2f%%\n\n",
				  $id, $seqlen, $gc_content;

	return $seqlen;

} # seq_stats

sub revcomplement() { ## Reverse Complement

	my ($seqref, $revseq);

	$seqref = shift;
	$revseq = shift;

	($$revseq = reverse $$seqref) =~ tr/acgtACGT/tgcaTGCA/;

} # revcomplement

sub init_codon_table() {

	%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

} # init_codon_table

sub translate() { ## Translate Sequence

	my ($seqref, $fr) = ($_[0], $_[1]); # @_

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

	my $seqlen = length($$seqref);

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

		my $protein_seq = '';

		for (my $i = $phase; $i < $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";

}; # translate


__END__
#
# REMEMBER: Anyting following the '__END__'
#           or '__DATA__' tokens is ignored
#
# 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
