BEGIN{
    $ENV{'PATH'}=$ENV{'PATH'}.":/local/workspace01/zhangju/script/perl";
    push(@INC,"/local/workspace01/zhangju/script/perl");
}

use strict;
use Bio::SeqIO;
use Bio::Annotation::Collection;
use Bio::Annotation::Comment;
 
my $input_file = "GenePRIMP.art";#"GenePRIMP.art";

my $root_dir="/local/workspace01/zhangju/wc1_annotation/curated";

my $embl_dir=$root_dir."/single_curated_embl";
my $gbk_dir=$root_dir."/single_curated_gbk";
my $new_embl=$root_dir."/single_new_embl";


my $seq_in  = Bio::SeqIO->new(
                              -format => 'genbank',
                              -file   => $input_file,
                              );
 
my $seq;

open(OUTFILE,">translation_issue.txt");

while( $seq = $seq_in->next_seq() ) {

    print OUTFILE "================\n";
    print OUTFILE "File - ".$seq->display_id."\n\n";
    my $embl_file=$embl_dir."/".$seq->display_id; #"/twc1a_scaffold00002.2";
    my $embl_in  = Bio::SeqIO->new(
                               -format => 'embl',
                               -file   => $embl_file,
			       );


    my $embl_seq=$embl_in->next_seq();
    for my $feat_object ($embl_seq->get_SeqFeatures) {          
       #print "primary tag: ", $feat_object->primary_tag, "\n"; 
       if ($feat_object->primary_tag eq "CDS"){
       	   my $start=$feat_object->location->start;
	   my $end=$feat_object->location->end;
	   my $complementary=$feat_object->strand;
	   my $seqobj=Bio::Seq->new( -seq => (substr $seq->seq(),$start-1,$end-$start+1),
	   			     -id => 'temp',
				     );
	   my $dnaobj=$seqobj;
	   if ($complementary==-1){
	   	$dnaobj=$seqobj->revcom;
	   }
	   my $protein=$dnaobj->translate(-codontable_id => 11);
	   my $pro_from_dna=$protein->seq();
				         
	   my $locus="";   
	   for my $tag ($feat_object->get_all_tags) { 
               #print $tag."\n";
               if ($tag eq "locus_tag"){

		   my @arr=$feat_object->get_tag_values($tag);
		   $locus=$arr[0]."\n";

	       }

               if ($tag eq "translation"){

		   my @arr1=$feat_object->get_tag_values($tag);
		   my $aa=$arr1[0];
		   if ($aa =~ m/\s/){
	               print OUTFILE $locus." tranlation with spaces, please verify\n";
		       #my @arr=$feat_object->get_tag_values($tag);
		       print OUTFILE $aa."\n";
		   }
		   if (not (($aa =~ m/\Q$pro_from_dna\E/) || ($pro_from_dna =~ m/\Q$aa\E/)) || abs(length($aa)-length($pro_from_dna))>5){
		       print OUTFILE $locus."\ttranslations mismatch, please verify\n\n";
		       print OUTFILE "stand - ".$complementary." coordinate (".$start."..".$end.")\n\n";
		       print OUTFILE $dnaobj->seq()."\n\n";
		       print OUTFILE "From dna direct translation: length - ".length($pro_from_dna)."  ".$pro_from_dna."\n\n";
		       print OUTFILE "Annotation translation: length - ".length($aa)."  ".$aa."\n\n";
		   }
	       }
	  }          
       }       
    }
}
