#!/usr/bin/perl
# snp_info outputs additional snp information is released under the
# terms of the GPL version 2, or any later version, at your option.
# See the file README and COPYING for more information.
# Copyright 2011 by Don Armstrong <don@donarmstrong.com>.


use warnings;
use strict;

BEGIN {$|=1};

use Getopt::Long;
use Pod::Usage;

=head1 NAME

snp_info - output snp information

=head1 SYNOPSIS

 [options]

 Options:
  --debug, -d debugging level (Default 0)
  --help, -h display this help
  --man, -m display manual

=head1 OPTIONS

=over

=item B<--nearest-gene>

Return the name of the nearest gene if the SNP isn't in a gene

=item B<--debug, -d>

Debug verbosity. (Default 0)

=item B<--help, -h>

Display brief usage information.

=item B<--man, -m>

Display this manual.

=back

=head1 EXAMPLES

echo 'rs1782034' | snp_info;

=cut


use vars qw($DEBUG);
use DBI;
use Data::Dumper;
use Scalar::Util qw(looks_like_number);


my %options = (debug           => 0,
	       help            => 0,
	       man             => 0,
	       verbose         => 0,
	       fix_pos         => 1,
	       nearest_gene    => 0,
	       build           => '135',
	       refseq          => '37_3',
	       );

GetOptions(\%options,
	   'verbose|v+',
	   'fix_pos|fix-pos!',
	   'nearest_gene|nearest-gene!',
	   'debug|d+','help|h|?','man|m');

pod2usage() if $options{help};
pod2usage({verbose=>2}) if $options{man};

$DEBUG = $options{debug};

my @USAGE_ERRORS;
# if (1) {
#      push @USAGE_ERRORS,"You must pass something";
# }

pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;


my $dbh = DBI->connect("dbi:Pg:service=snp",'','',{AutoCommit => 0}) or
    die "Unable to connect to database: ".$DBI::errstr;

my $real_snp_sth = $dbh->prepare(<<'END') // die "Unable to prepare realsnp statement: ".$dbh->errstr;
SELECT DISTINCT rscurrent FROM rsmergearch WHERE rslow=? OR rshigh=?
END

my $sth_snp_history = $dbh->prepare(<<'END') // die "Unable to prepare snp_history statement: ".$dbh->errstr;
SELECT * FROM snphistory WHERE snp_id=?;
END

my $snp_chr_pos_sth = $dbh->prepare(<<"END") // die "Unable to prepare snpinfo statement: ".$dbh->errstr;
SELECT scpr.snp_id AS snp_id,
       scpr.chr AS chr,
       scpr.pos AS pos,
       scl.phys_pos_from AS pos2,
       scpr.orien AS orien,
       scl.allele AS ref,
       uv.var_str AS var_str,
       ruv.var_str AS rev_var_str,
       ml.locus_id AS uid,
       ml.locus_symbol AS symbol,
       gitn.gene_name AS description,
       scl.aln_quality AS alignment_quality
       FROM snp s
         JOIN b$options{build}_snpchrposonref_$options{refseq} scpr ON s.snp_id=scpr.snp_id
         JOIN b$options{build}_snpcontigloc_$options{refseq} scl ON scpr.snp_id=scl.snp_id
         JOIN b$options{build}_contiginfo_$options{refseq} ci ON scl.ctg_id = ci.ctg_id
         LEFT OUTER JOIN b$options{build}_snpcontiglocusid_$options{refseq} ml ON s.snp_id=ml.snp_id
         LEFT OUTER JOIN geneidtoname gitn ON ml.locus_id=gitn.gene_id
         JOIN univariation uv ON s.univar_id=uv.univar_id
         JOIN univariation ruv ON uv.rev_univar_id=ruv.univar_id
WHERE (ci.group_term LIKE 'GRCh%' OR ci.group_term='Primary_Assembly') AND s.snp_id=? ORDER BY scl.aln_quality;
END

my $sth = $dbh->prepare(<<"END") // die "Unable to prepare snpinfo statement: ".$dbh->errstr;
SELECT snp_id FROM b$options{build}_snpchrposonref_$options{refseq} WHERE chr=? AND pos=?;
END

my $sth_mrna_ss = $dbh->prepare(<<'END') // die "Unable to prepare mrna_ss statement: ".$dbh->errstr;
SELECT * FROM mrna WHERE start<=? AND stop>=? AND chr=? ORDER BY variant;
END

my $sth_cds_ss = $dbh->prepare(<<'END') // die "Unable to prepare cds_ss statement: ".$dbh->errstr;
SELECT * FROM cds WHERE start<=? AND stop>=? AND chr=? ORDER BY variant;
END

my $sth_mrna_cds = $dbh->prepare(<<'END') // die "Unable to prepare mrna_cds statement: ".$dbh->errstr;
SELECT * FROM mrna_cds_table WHERE cds_gi=?;
END

my $sth_allele_freq = $dbh->prepare(<<'END') // die "Unable to prepare allele_freq statement: ".$dbh->errstr;
SELECT ssl.snp_id AS snp_id,
       si.subsnp_id AS subsnp_id,
       si.chr AS chr,
       si.submitted_ind_id AS submitted_ind_id,
       ga.chr_num AS chr_num,
       ga.rev_flag AS rev_flag,
       a.allele AS allele
  FROM snpsubsnplink ssl
    JOIN subind si on ssl.subsnp_id=si.subsnp_id
    JOIN gtyallele ga on si.gty_id=ga.gty_id
    JOIN allele a on ga.fwd_allele_id = a.allele_id
    WHERE ssl.snp_id=? AND
          si.chr=? AND ga.rev_flag = '0' AND
          a.allele <> 'N';
END

my $sth_rsid_info = $dbh->prepare(<<'END') // die "Unable to prepare rsid info statement: ".$dbh->errstr;
SELECT * FROM cds WHERE start<=? AND stop>=? AND chr=? ORDER BY variant;
END

my $sth_near_gene_start = $dbh->prepare(<<'END') // die "Unable to prepare near gene start statement: ".$dbh->errstr;
SELECT cds.start AS start,
       cds.stop AS stop,
       cds.gene AS symbol,
       gitn.gene_name AS description
    FROM cds
      LEFT OUTER JOIN geneidtoname gitn ON cds.geneid=gitn.gene_id
    WHERE cds.chr = ? AND cds.start > ?
    ORDER BY cds.start ASC LIMIT 1;
END

my $sth_near_gene_stop = $dbh->prepare(<<'END') // die "Unable to prepare near gene stop statement: ".$dbh->errstr;
SELECT cds.start AS start,
       cds.stop AS stop,
       cds.gene AS symbol,
       gitn.gene_name AS description
    FROM cds
      LEFT OUTER JOIN geneidtoname gitn ON cds.geneid=gitn.gene_id
    WHERE cds.chr = ? AND cds.stop < ?
    ORDER BY cds.stop DESC LIMIT 1;
END

my $prot_table;

while (<DATA>) {
    chomp;
    my ($codon,$aa,$aa_med,$aa_long) = split /\t/;
    $prot_table->{$codon} = $aa;
}

my @headers = qw(id chr pos ref alt);
my @row_keys_to_initialize = qw(gene nonsyn aa_ref aa_alt codon_ref codon_alt cdspos protpos loc num var acc phase reffreq altfreq refaltn);
print join("\t",@headers,'orig_id',@row_keys_to_initialize)."\n";
while (<>) {
    print $_ and next if /^#/;
    chomp;
    my @row = split /\t/, $_;
    my %row;
    @row{@headers} = @row[0..$#headers];
    $row{orig_id} = $row{id};
    # this is taken directly from vcf_fill_rsid; we want it for the time being anyway
    if ($row{id} !~ /^(?:rs|ss)?\d+(?:,\s*(?:rs|ss)?\d+)*$/) {
	my $new_id = find_rsid(chr => fixup_chr($row{chr}),
			       # snps are 0 indexed
			       pos => $row{pos} - 1,
			       dbh => $dbh,
			       sth => $sth,
			      );
	if ($new_id !~ /^(?:[\.\?])$/) {
	    $row{id} = "rs".$new_id;
	}
    }
    if ($row{id} =~ /(rs|ss)(\d+)/) {
	my $type = $1;
	my $id = $2;
	$id = find_real_snp_id(id => $id,
			       dbh => $dbh,
			       sth => $real_snp_sth,
			      );
	if (not defined $id and
	    is_retired_snp(id  => $id,
			   dbh => $dbh,
			   sth => $sth_snp_history,
			  )) {
	    if (defined $row{chr} and defined $row{pos}) {
		$row{id} = "$row{chr}-$row{pos}";
	    }
	}
	else {
	    $row{id} = "$type$id";
	}
    }

    if (defined $row{id} and ($options{fix_pos} or
			      not defined $row{pos} or
			      not defined $row{chr} or
			      not defined $row{ref} or
			      not defined $row{alt})) {
	my %info = find_chr_pos(id => [$row{id} =~ /(?:rs|ss)(\d+)/]->[0],
				dbh => $dbh,
				sth => $snp_chr_pos_sth,
			       );
	for my $inf (qw(pos chr ref alt)) {
	    $row{$inf} = $info{$inf} if $options{fix_pos} and defined $info{$inf} or not defined $row{$inf};
	}
	# handle partially aligned snps
	if ($options{fix_pos} and not defined $row{pos}) {
	    $row{pos} = $info{pos2};
	}
	if (not defined $row{pos} or not defined $row{chr}) {
	    if (defined $row{pos}) {
		print STDERR "$row{id} is ambiguously placed\n";
	    }
	    else {
		print STDERR "Unable to find a position/chr for $row{id}\n";
	    }
	}
    }
    my $mrna_information = mrna_information(dbh => $dbh,
					    sth => $sth_mrna_ss,
					    # these are all 1 indexed
					    pos => $row{pos},
					    chr => fixup_chr($row{chr}),
					    ref => $row{ref},
					    alt => $row{alt},
					    cds_sth => $sth_cds_ss,
					    mrna_cds_sth => $sth_mrna_cds,
					   );
    @row{@row_keys_to_initialize} = ('') x @row_keys_to_initialize;
    if (keys %{$mrna_information}) {
	# strip out existing information in $row{info}
	my @extra_info;
	$row{gene}     = $mrna_information->{gene};
	$row{nonsyn}   = 1 if defined $mrna_information->{nonsyn} and $mrna_information->{nonsyn};
	$row{aa_ref}    = $mrna_information->{orig_aa} if defined $mrna_information->{orig_aa};
	$row{aa_alt}    = $mrna_information->{new_aa} if defined $mrna_information->{new_aa};
	$row{codon_ref} = $mrna_information->{orig_codon} if defined $mrna_information->{orig_codon};
	$row{codon_alt} = $mrna_information->{new_codon} if defined $mrna_information->{new_codon};
	$row{cdspos}   = $mrna_information->{cds_pos} if defined $mrna_information->{cds_pos};
	$row{protpos}  = $mrna_information->{prot_pos} if defined $mrna_information->{prot_pos};
	$row{loc}      = $mrna_information->{exon}?'exon':'intron' if defined $mrna_information->{gene};
	$row{num}      = $mrna_information->{number};
	$row{var}      = $mrna_information->{variant} if defined $mrna_information->{variant};
	$row{acc}      = $mrna_information->{accession};
	$row{phase}    = $mrna_information->{phase} if defined $mrna_information->{phase};
    }
    if ((not defined $row{gene} or not length $row{gene})
	and $options{nearest_gene}) {
	$row{gene} = find_nearest_gene(dbh => $dbh,
				       start_sth => $sth_near_gene_start,
				       stop_sth => $sth_near_gene_stop,
				       chr => fixup_chr($row{chr}),
				       pos => $row{pos},
				      );
    }
    my $allele_freq_information = allele_freq_information(dbh => $dbh,
							  sth => $sth_allele_freq,
							  id  => [$row{id} =~ /(?:rs|ss)(\d+)/]->[0],
							  # these are all 1 indexed
							  pos => $row{pos},
							  chr => fixup_chr($row{chr}),
							  ref => $row{ref},
							  alt => $row{alt},
							 ) if defined $row{id} and $row{id} =~ /^(?:rs|ss)/;
    if (keys %{$allele_freq_information}) {
	# strip out existing information in $row{info}
	if (exists $allele_freq_information->{alleles} and keys %{$allele_freq_information->{alleles}}) {
	    $row{reffreq} = sprintf('%.3f',allele_freq($allele_freq_information->{alleles},$row{ref}));
	    $row{altfreq} = sprintf('%.3f',allele_freq($allele_freq_information->{alleles},$row{alt}));
	    $row{refaltn} = $allele_freq_information->{alleles}{total};
	}
    }
    print join("\t",map {defined $_?$_:'NA'} @row{@headers,'orig_id',@row_keys_to_initialize},$#row > $#headers?@row[@headers..$#row]:())."\n";
    if ($options{verbose} > 0 and $. % 100 == 0) {
	print STDERR "At line $.\n";
    }
}

sub find_real_snp_id {
    my %param = @_;

    my %info;
    my $rv = $param{sth}->execute($param{id},$param{id}) //
	die "Unable to execute statement properly: ".$param{dbh}->errstr;
    my ($real_id) = map {ref $_ ?@{$_}:()} map {ref $_ ?@{$_}:()} $param{sth}->fetchall_arrayref([0]);
    if ($param{sth}->err) {
	print STDERR $param{sth}->errstr;
	$param{sth}->finish;
	return $param{id};
    }
    $param{sth}->finish;
    return $real_id // $param{id};
}

sub find_chr_pos {
    my %param = @_;

    my %info;
    my $rv = $param{sth}->execute($param{id}) //
	die "Unable to execute statement properly: ".$param{dbh}->errstr;
    my $results = $param{sth}->fetchall_arrayref({});
    if ($param{sth}->err) {
	print STDERR $param{sth}->errstr;
	$param{sth}->finish;
	return %info;
    }
    for my $result (@{$results}) {
	%info = %{$result};
    }
    if (keys %info) {
	if ($info{chr} eq 'NotOn') {
	    undef $info{chr};
	}
	if ($info{orien}) {
	    $info{orien_var_str} = $info{rev_var_str};
	}
	else {
	    $info{orien_var_str} = $info{var_str};
	}
	my @alleles = split qr{/},$info{orien_var_str};
	($info{alt}) = grep {$_ ne $info{ref}} @alleles;
	# fix $row{pos}
	$info{pos}+=1;
    }
    $param{sth}->finish;
    return %info;
}

sub find_nearest_gene {
    my %param = @_;
    return undef unless defined $param{chr} and defined $param{pos};

    my %possible_genes;

    for my $st (qw(start stop)) {
	my $rv = $param{"${st}_sth"}->execute($param{chr},$param{pos}) //
	    die "Unable to execute nearest gene $st statement properly: ".$param{dbh}->errstr;
	my $results = $param{"${st}_sth"}->fetchall_arrayref({});
	if ($param{"${st}_sth"}->err) {
	    print STDERR $param{"${st}_sth"}->errstr;
	    $param{"${st}_sth"}->finish;
	    return '';
	}
	$possible_genes{$st} = $results->[0]
    }
    for my $st (qw(start stop)) {
	# opposite of start/stop (stop/start)
	my $st_op = $st eq q(start) ? q(stop) : q(start);
	if (not defined $possible_genes{$st}{start}) {
	    if (defined $possible_genes{$st_op}{start}) {
		return $possible_genes{$st_op}{symbol}
	    }
	    else {
		return undef;
	    }
	}
    }
    my $gene;
    # if the start of the start gene is closer to the position of the
    # snp than the stop of the stop gene, it's the nearest gene;
    # otherwise, the stop gene is
    if ($possible_genes{start}{start} - $param{pos} <
	$param{pos} - $possible_genes{stop}{stop}
       ) {
	return $possible_genes{start}{symbol};
    }
    else {
	return $possible_genes{stop}{symbol};
    }
}


sub mrna_information{
    my %param = @_;

    my $rv = $param{sth}->execute($param{pos},$param{pos},$param{chr}) or
	die "Unable to execute statement properly: ".$param{dbh}->errstr;

    my $results = $param{sth}->fetchall_arrayref({});
    if ($param{sth}->err) {
	print STDERR $param{sth}->errstr;
	return {};
    }
    my $chosen_result;
    for my $result (@{$results}) {
	if (not defined $chosen_result) {
	    $chosen_result = $result;
	    next;
	}
	if (not defined $chosen_result->{variant} or
	    (defined $result->{variant} and
	     $chosen_result->{variant} > $result->{variant})) {
	    $chosen_result = $result;
	    next;
	}
    }
    if (defined $chosen_result and $chosen_result->{exon}) {
	# try to figure out the cds if there is mrna information and
	# we're in an exon
	my $cds_rv = $param{cds_sth}->execute($param{pos},$param{pos},$param{chr}) or
	    die "Unable to execute statement properly: ".$param{dbh}->errstr;
	my $cds_results = $param{cds_sth}->fetchall_arrayref({});
	if ($param{cds_sth}->err) {
	    print STDERR $param{cds_sth}->errstr;
	}
	else {
	    print STDERR Data::Dumper->Dump([$cds_results],[qw(cds_results)]) if $DEBUG;
	    my $chosen_cds;
	    for my $cds_result (@{$cds_results}) {
		if (not defined $chosen_cds) {
		    $chosen_cds = $cds_result;
		    next;
		}
		if (not defined $cds_result->{variant} or
		    (defined $chosen_cds->{variant} and
		     $chosen_cds->{variant} ge $cds_result->{variant})) {
		    $chosen_cds = $cds_result;
		    next;
		}
	    }
	    if (defined $chosen_cds and keys %{$chosen_cds}) {
		my $comp_offset = $chosen_cds->{complement}?-1:1;
		$chosen_result->{in_cds} = 1;
		$chosen_result->{cds_pos} =
		    ($chosen_cds->{complement}?$chosen_cds->{stop}-$param{pos}:$param{pos}-$chosen_cds->{start}) +
			$chosen_cds->{position};
		$chosen_result->{gene} = $chosen_cds->{gene};
		$chosen_result->{pos} = $chosen_cds->{pos};
		$chosen_result->{number} = $chosen_cds->{number};
		$chosen_result->{variant} = $chosen_cds->{variant};
		$chosen_result->{acc} = $chosen_cds->{accession};
		$chosen_result->{prot_pos} = sprintf("%.1f",$chosen_result->{cds_pos}/ 3);
		$chosen_result->{phase} = ($chosen_cds->{phase} + $param{pos} - $chosen_cds->{start}) % 3;
		print STDERR Data::Dumper->Dump([$chosen_cds,$chosen_result],[qw(chosen_cds chosen_result)]) if $DEBUG;
		my $mrna_cds_rv = $param{mrna_cds_sth}->execute($chosen_cds->{gi}) or
		    die "Unable to execute statement properly: ".$param{dbh}->errstr;
		my $mrna_cds_result = $param{mrna_cds_sth}->fetchrow_hashref();
		if (defined $mrna_cds_result and keys %{$mrna_cds_result}) {
 		    print STDERR Data::Dumper->Dump([$mrna_cds_result,$chosen_result,\%param],
						    [qw(mrna_cds_result chosen_result *param)]) if $DEBUG;
		    my $codon = substr $mrna_cds_result->{sequence},$chosen_result->{cds_pos} - (($chosen_result->{cds_pos}) % 3)  ,3;
		    my $orig_aa = convert_to_prot_seq($codon);
		    my $modified_codon = $codon;
		    substr($modified_codon,$chosen_result->{cds_pos} % 3,1)=
			($chosen_cds->{complement}?dna_complement($param{alt}):$param{alt});
		    my $new_aa = convert_to_prot_seq($modified_codon);
  		    print STDERR Dumper($mrna_cds_result,$chosen_result,$chosen_cds,{codon_pos => $chosen_result->{cds_pos} - ($chosen_result->{cds_pos} % 3),
 										     codon=>$codon,modified_codon=>$modified_codon,pos=>$param{pos}}) if $DEBUG;
		    if ($DEBUG > 1) {
			my $a = -1;
			print STDERR map {$a++; ($a*3)." $_\n"} $mrna_cds_result->{sequence} =~ m/.../g;
		    }
 		    $chosen_result->{orig_aa} = $orig_aa;
		    $chosen_result->{new_aa} = $new_aa;
		    $chosen_result->{new_codon} = $modified_codon;
		    $chosen_result->{orig_codon} = $codon;
		    $chosen_result->{nonsyn} = $new_aa ne $orig_aa ? 1 : 0;
		}
	    }
	}
    }
    return $chosen_result // {};
}


sub fixup_chr {
    my ($chr) = @_;
    return $chr if not defined $chr;

    $chr =~ s/^NC_00000?//;
    $chr =~ s/^chr//;
    if (looks_like_number($chr)) {
	if ($chr == 23) {
	    $chr = "X"
	}
	elsif ($chr == 24) {
	    $chr = "Y"
	}
    }
    return ($chr);
}

sub allele_freq_information{
    my %param = @_;

    return {} if $param{id} !~ /^\d+$/;

    my $rv = $param{sth}->execute($param{id},$param{chr}) or
	die "Unable to execute statement properly: ".$param{dbh}->errstr." id:".$param{id}." chr:".$param{chr};

    my $results = $param{sth}->fetchall_arrayref({});
    if ($param{sth}->err) {
	print STDERR $param{sth}->errstr;
	return {};
    }
    my %pop_alleles;
    my %alleles;
    my %seen_indiv;
    my %indivs;
    for my $result (@{$results}) {
	next if exists $seen_indiv{$result->{submitted_ind_id}.'.'.$result->{chr_num}};
	$indivs{$result->{submitted_ind_id}}{'chr_'.$result->{chr_num}} = $result->{allele};
	if (defined $result->{loc_pop_id}) {
	    $indivs{$result->{submitted_ind_id}}{loc_pop_id} = $result->{loc_pop_id};
	    $indivs{$result->{submitted_ind_id}}{pop_handle} = $result->{pop_handle};
	    $pop_alleles{$result->{loc_pop_id}}{$result->{allele}}++;
	    $pop_alleles{$result->{loc_pop_id}}{total}++;
	}
	$alleles{$result->{allele}}++;
	$alleles{total}++;
	$seen_indiv{$result->{submitted_ind_id}.'.'.$result->{chr_num}} = 1;
    }
    my %genotypes;
    my %pop_genotypes;
    for my $indiv (values %indivs) {
	if (defined $indiv->{chr_1} and defined $indiv->{chr_2}) {
	    my $genotype = join('/',
				sort($indiv->{chr_1},
				     $indiv->{chr_2}));
	    if (defined $indiv->{loc_pop_id}) {
		$pop_genotypes{$indiv->{loc_pop_id}}{$genotype}++;
		$pop_genotypes{$indiv->{loc_pop_id}}{total}++;
	    }
	    $genotypes{$genotype}++;
	    $genotypes{total}++;
	}
    }
    return {pop_alleles   => \%pop_alleles,
	    pop_genotypes => \%pop_genotypes,
	    genotypes     => \%genotypes,
	    alleles       => \%alleles,
	   };
}

sub allele_freq{
    my ($freq_table,$allele) = @_;
    my $flipped_allele = $allele;
    $flipped_allele =~ tr/ACTG/TGAC/;
    if (defined $allele) {
	if (exists $freq_table->{$allele} and defined $freq_table->{$allele}) {
	    return $freq_table->{$allele}/$freq_table->{total};
	}
	elsif (exists $freq_table->{$flipped_allele} and defined $freq_table->{$flipped_allele}) {
	    return $freq_table->{$flipped_allele}/$freq_table->{total};
	}
	return 0;
    }
    my @alleles = grep {$_ ne 'total'} keys %{$freq_table};
    return join(',',map {$_.':'.(defined $freq_table->{$_} ?
				 $freq_table->{$_}/$freq_table->{total} : 0)} @alleles);
}


sub find_rsid {
    my %param = @_;

    my $rv = $param{sth}->execute($param{chr},$param{pos}) //
	die "Unable to execute statement properly: ".$param{dbh}->errstr;
    my @snp_ids = grep {$_ =~ /^\d+$/}
	map {ref $_ eq 'ARRAY'? @{$_}:$_}
	    map {ref $_ eq 'ARRAY'? @{$_}:$_}
		$sth->fetchall_arrayref([0]);
    if ($sth->err) {
	print STDERR $sth->errstr;
	return '?';
    }
    if (not @snp_ids) {
	return '.';
    }
    return join(',',@snp_ids);
}

sub dna_complement{
    my ($seq) = @_;
    return join('',map {tr /ATCG/TAGC/; $_} reverse split //,$seq);
}

sub convert_to_prot_seq{
    my ($mrna) = @_;
    $mrna =~ s/U/T/g;
    my $prot = '';
    for my $pos (0..(length($mrna) / 3 - 1)) {
	my $codon = substr($mrna,$pos*3,3);
	if (not exists $prot_table->{$codon}) {
	    # this isn't quite right, as we should really be handling
	    # ambiguities which do not result in different amino
	    # acids, but for the time being, do it this way.
	    $prot .= 'N';
	}
	else {
	    $prot .= $prot_table->{$codon};
	}
    }
    return $prot;
}


__DATA__
ATT	I	Ile	Isoleucine
ATC	I	Ile	Isoleucine
ATA	I	Ile	Isoleucine
CTT	L	Leu	Leucine
CTC	L	Leu	Leucine
CTA	L	Leu	Leucine
CTG	L	Leu	Leucine
TTA	L	Leu	Leucine
TTG	L	Leu	Leucine
GTT	V	Val	Valine
GTC	V	Val	Valine
GTA	V	Val	Valine
GTG	V	Val	Valine
TTT	F	Phe	Phenylalanine
TTC	F	Phe	Phenylalanine
ATG	M	Met	Methionine
TGT	C	Cys	Cysteine
TGC	C	Cys	Cysteine
GCT	A	Ala	Alanine
GCC	A	Ala	Alanine
GCA	A	Ala	Alanine
GCG	A	Ala	Alanine
GGT	G	Gly	Glycine
GGC	G	Gly	Glycine
GGA	G	Gly	Glycine
GGG	G	Gly	Glycine
CCT	P	Pro	Proline
CCC	P	Pro	Proline
CCA	P	Pro	Proline
CCG	P	Pro	Proline
ACT	T	Thr	Threonine
ACC	T	Thr	Threonine
ACA	T	Thr	Threonine
ACG	T	Thr	Threonine
TCT	S	Ser	Serine
TCC	S	Ser	Serine
TCA	S	Ser	Serine
TCG	S	Ser	Serine
AGT	S	Ser	Serine
AGC	S	Ser	Serine
TAT	Y	Tyr	Tyrosine
TAC	Y	Tyr	Tyrosine
TGG	W	Trp	Tryptophan
CAA	Q	Gln	Glutamine
CAG	Q	Gln	Glutamine
AAT	N	Asn	Asparagine
AAC	N	Asn	Asparagine
CAT	H	His	Histidine
CAC	H	His	Histidine
GAA	E	Glu	Glutamic acid
GAG	E	Glu	Glutamic acid
GAT	D	Asp	Aspartic acid
GAC	D	Asp	Aspartic acid
AAA	K	Lys	Lysine
AAG	K	Lys	Lysine
CGT	R	Arg	Arginine
CGC	R	Arg	Arginine
CGA	R	Arg	Arginine
CGG	R	Arg	Arginine
AGA	R	Arg	Arginine
AGG	R	Arg	Arginine
TAA	*	*	Stop codon
TAG	*	*	Stop codon
TGA	*	*	Stop codon