#!/usr/bin/perl

# this program is used to translate the cDNA sequence into amino acids using frame +1.
#
####################################################################
# Program name: Translation+1.1.pl
# Version 1.1 Revised
# written by Jack Min @copy right, on August 9, 2003

#######################

########################################################################################
# In the output amino acid sequences, X: unknow amino acids from translating nucleotide such as "N" or "n";
# the sign '*' is the stop codon translation.
#
# the amino acid translation start from the frame value for 5'-strand.
#
################ Body of the code #######################################

# the number of argument must be 2

if ($#ARGV != 1)
{
#	print "\nPlease check you input. The number of arguments is not correct!\n\n";
	print "Usage: Translation+.pl Seq-Fasta-file ProteinFile\n\n";
	print "Input:\n";
	print "		Seq-Fasta - DNA sequences in FASTA format;\n";
	print "Output:\n";
	print "		ProteinOutput - output of translated protein sequences\n\n";
	exit;
}

$query = "" ;  #query sequece: contig number or clone number

$frame =""; #used for translation

$head = 0; #false
$hit = 0;#false

@proteinArray = (); # to store translated sequence
##################################################
$result = "";
$proteinName = "";

$dnaFile = $ARGV[0];
$proteinSeqFile = $ARGV[1];

print "In progress ...\n";
$midLine = 0; #the line below the "Query:" to show identity and similarity
###################################################

$position1 = 0; # blank space
$position2 = 0; # delimiter position including ", (, ["
$noHitCount = 0;
$queryCount = 0;
$mass = 0;
$signal = '';

%seqH = (); # a hash for storing seqID (key) and seq (value)
$line =();
open (INPUT, $dnaFile) || die ("Cannot open file $dnaFile\n\n");

$seqFlag =0; #to indicate the line is header or seqs, used for seq in multiple lines
$value ='';
$key = '';
$count =0; #to indicate the first ">"
@idArray = ();
while (<INPUT>)
{
	$line = $_; #$_, default input variable
	chomp ($line);
	if ($line =~ /^>/)
	{
		$count++;
		if ($count >1)
		{
			$seqH{$key} = $value;
			# print $key."\n".$value."\n";

			$value ='';
			$key ='';
		}
		$seqFlag = 1;
		@tmpA = split(/ /, $line);
		$key = substr($tmpA[0], 1);
		push (@idArray, $key."\n");
	}
	else
	{
		if ($seqFlag ==1)
		{
			$value .= $line;
		}
	}
}
$seqH{$key} = $value;
close (INPUT);

foreach $entry (@idArray)
{
	chomp ($entry);
	#frame +1
	$frame = 1;
	$sequence = $seqH{$entry};
	$protein = dna2protein($frame, $sequence); #to obtain protein seq
	$mass = int (getMass ($protein));
	push (@proteinArray, ">".$entry."\t+1\t".$mass."\n".$protein."\n");
	# add code for frame +2, frame +3, frame -1, -2, -3. You can use the same subroutine for - frames, but need to 
	# use reverse-complementary strand	

#push (@dbArray, $entry."\t+1"."\t".$mass."\t".$protein."\n");
}

open (OUT, ">$proteinSeqFile");
print OUT (@proteinArray);
close (OUT);

$f2 = $proteinSeqFile.".db";
open (OUT, ">$f2");
print OUT (@dbArray);
close (OUT);

print "The work is done!!\n";

exit;
#########################################
# A subroutine to translate a DNA 3 character codon to an amino acid
#
# The code is adapted by Jack Min from Beginning Perl for Bioinformatics, Authour James Tisdall, Page 159
###############################################################################
sub codon2aa
{
	use strict;
	my ($codon) = @_;

        if ( $codon =~ /GC./i)        { return 'A' }    # Alanine
     elsif ( $codon =~ /TG[TC]/i)     { return 'C' }    # Cysteine
     elsif ( $codon =~ /GA[TC]/i)     { return 'D' }    # Aspartic Acid
     elsif ( $codon =~ /GA[AG]/i)     { return 'E' }    # Glutamic Acid
     elsif ( $codon =~ /TT[TC]/i)     { return 'F' }    # Phenylalanine
     elsif ( $codon =~ /GG./i)        { return 'G' }    # Glycine
     elsif ( $codon =~ /CA[TC]/i)     { return 'H' }    # Histidine
     elsif ( $codon =~ /AT[TCA]/i)    { return 'I' }    # Isoleucine
     elsif ( $codon =~ /AA[AG]/i)     { return 'K' }    # Lysine
     elsif ( $codon =~ /TT[AG]|CT./i) { return 'L' }    # Leucine
     elsif ( $codon =~ /ATG/i)        { return 'M' }    # Methionine
     elsif ( $codon =~ /AA[TC]/i)     { return 'N' }    # Asparagine
     elsif ( $codon =~ /CC./i)        { return 'P' }    # Proline
     elsif ( $codon =~ /CA[AG]/i)     { return 'Q' }    # Glutamine
     elsif ( $codon =~ /CG.|AG[AG]/i) { return 'R' }    # Arginine
     elsif ( $codon =~ /TC.|AG[TC]/i) { return 'S' }    # Serine
     elsif ( $codon =~ /AC./i)        { return 'T' }    # Threonine
     elsif ( $codon =~ /GT./i)        { return 'V' }    # Valine
     elsif ( $codon =~ /TGG/i)        { return 'W' }    # Tryptophan
     elsif ( $codon =~ /TA[TC]/i)     { return 'Y' }    # Tyrosine
     elsif ( $codon =~ /TA[AG]|TGA/i) { return '*' }    # Stop
     else { return 'X'} # for nucleotide as 'N' or n or others not 'AGCT'

 }


############################################
# A subroutine to translate DNA sequence int a protein
#
# the code is adapted by Jack Min from the same source as above, on Page 165
###############################################
sub dna2protein
{
	my ($frame, $dna) = @_;
	use strict;
	my ($protein) = '';
	for (my $i =$frame-1; $i < (length($dna) -2); $i +=3)
	{
		$protein .=codon2aa (substr($dna, $i, 3));
	}

	return $protein;
}

 #############################
 # A subroutine to match molecular weight of each amino acid residue
 #
 #############################
sub AAMass
{
	use strict;
	my ($AA) = @_;
	if ($AA eq 'A') {return 71.1}
	elsif ($AA eq 'R') {return 156.2}
	elsif ($AA eq 'N') {return 114.1}
	elsif ($AA eq 'D') {return 115.1}
	elsif ($AA eq 'C') {return 103.1}
	elsif ($AA eq 'Q') {return 128.1}
	elsif ($AA eq 'E') {return 129.1}
	elsif ($AA eq 'G') {return 57.1}
	elsif ($AA eq 'H') {return 137.2}
	elsif ($AA eq 'I') {return 113.2}
	elsif ($AA eq 'L') {return 113.2}
	elsif ($AA eq 'K') {return 128.2}
	elsif ($AA eq 'M') {return 131.2}
	elsif ($AA eq 'F') {return 147.2}
	elsif ($AA eq 'P') {return 97.1}
	elsif ($AA eq 'S') {return 87.1}
	elsif ($AA eq 'T') {return 101.1}
	elsif ($AA eq 'W') {return 186.2}
	elsif ($AA eq 'Y') {return 163.2}
	elsif ($AA eq 'V') {return 99.1}
	elsif ($AA eq 'X') {return 0.0}
	else {return 0.0}
}

#################################################
#
# A subroutine to get the mass of a predicted peptide
#
##################################################
sub getMass
{
	use strict;
	my ($peptide) = @_;
	my ($mass) = 0;
	my ($aa) = '';
	my ($stopFlag) = 0;
	my ($stopAfterStartFlag) = 0;
	my ($startFlag) =0;
	my ($length) = 0;
	$length = length($peptide);



	for (my ($i) = 0; $i < (length($peptide)); $i++)
	{
		$aa = substr($peptide, $i, 1);
		$mass += AAMass($aa);
	}

	return $mass;
}

