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

use Bio::SeqIO;
use Bio::Perl;
use Bio::Seq;
use Bio::Factory::FTLocationFactory;
use Bio::DB::GenPept;
use Bio::DB::GenBank;
use List::Util qw[min max];
use Bio::SeqUtils;
use Bio::SeqFeature::Generic;
use List::MoreUtils qw(uniq);

sub extract_CSterc{
	my %gene=();
	open transF,">".$gbk_file."_transcripts.txt";
	open genomeF,">".$gbk_file."_genome.fasta";
	while(my $seq_object=$inseq->next_seq){
		#print "primary tag: ", $feat_object->primary_tag, "\n"; 
		print "Pause";
		print genomeF ">".$seq_object->id."\n".$seq_object->seq."\n";
		my $contig_id=$seq_object->id;
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			
			my $start=0;
			my $end=0;
			if ($feat_object->has_tag("locus_tag")){
				my $pos1=$feat_object->start;
				my $pos2=$feat_object->end;
				$start=min($pos1,$pos2);
				$end=max($pos1,$pos2);
				my $key=$contig_id.":".($start-1)."-".$end;
				my $tag=($feat_object->get_tag_values("locus_tag"))[0];
				my $function="N/A";
				if($feat_object->has_tag("product")){
					$function=($feat_object->get_tag_values("product"))[0];
				}elsif($feat_object->has_tag("function")){
					$function=($feat_object->get_tag_values("function"))[0];
				}elsif($feat_object->has_tag("note")){
					$function=($feat_object->get_tag_values("note"))[0];
				}elsif($feat_object->has_tag("pseudogene")){
					$function="pseudogene"
				}
				$gene{$tag}=$feat_object->primary_tag."\t".$tag."\t".length($feat_object->spliced_seq->seq)."\t".$function;
			}
		}
	}
	
	for(keys %gene){
		print transF $gene{$_}."\n";
	}
	close(transF);
	genomeF.close();
}

sub extract_gd1_genePRIMP{
	open FPKMLocF, "<fpkm_locus.txt";
	my %gene=();
	foreach $line(<FPKMLocF>){
    	chomp($line);
    	$gene{$line}="";
    
	}
	open transF,">".$gbk_file."_transcripts.txt";
	
	
	while(my $seq_object=$inseq->next_seq){
		#print "primary tag: ", $feat_object->primary_tag, "\n"; 
		print "Pause";
		my $contig_id=$seq_object->id;
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			
			my $start=0;
			my $end=0;
			if ($feat_object->has_tag("locus_tag")){
				my $pos1=$feat_object->start;
				my $pos2=$feat_object->end;
				$start=min($pos1,$pos2);
				$end=max($pos1,$pos2);
				my $key=$contig_id.":".($start-1)."-".$end;
				my $tag=($feat_object->get_tag_values("locus_tag"))[0];
				my $function="N/A";
				if($feat_object->has_tag("product")){
					$function=($feat_object->get_tag_values("product"))[0];
				}elsif($feat_object->has_tag("function")){
					$function=($feat_object->get_tag_values("function"))[0];
				}elsif($feat_object->has_tag("note")){
					$function=($feat_object->get_tag_values("note"))[0];
				}elsif($feat_object->has_tag("pseudogene")){
					$function="pseudogene"
				}
				if(exists $gene{$key}){
					$gene{$key}=$feat_object->primary_tag."\t".$tag."\t".length($feat_object->spliced_seq->seq)."\t".$function;
				}
			}
		}
	}

	for(keys %gene){
		print transF $_."\t".$gene{$_}."\n";
	}
	close(transF);
}

sub extract_genes_canola{

	open transF,">".$gbk_file."_gene.fna";
	
	my %gene=();
	my @product=();
	my $gene_tag="tag";
	my $gene_seq="";
	while(my $seq_object=$inseq->next_seq){
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if ($feat_object->primary_tag eq "gene"){
				if($gene_tag ne "tag"){
					print transF ">".$gene_tag."|".join(";",uniq(@product))."\n".$gene_seq."\n"
				}
				@product=();
				$gene_tag=($feat_object->get_tag_values("locus_tag"))[0];
				$gene_seq=$feat_object->spliced_seq->seq;
			}
			if($feat_object->has_tag("product")){
	        	$product1=($feat_object->get_tag_values("product"))[0];
	        	push(@product, join('_',split(" ", $product1)));
			}
		}
		print transF ">".$gene_tag."|".join(";",@product)."\n".$gene_seq."\n"
		
	}

	close(transF);
}

sub extract_genes{

	open transF,">".$gbk_file."_gene.fna";
	
	my %gene=();
	
	while(my $seq_object=$inseq->next_seq){
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			
			if ($feat_object->primary_tag eq "CDS"){
				my $tag=($feat_object->get_tag_values("locus_tag"))[0];
				my $gene_name="";
				my $product="";
				if($feat_object->has_tag("product")){
		        	$product=($feat_object->get_tag_values("product"))[0];
				}
				my $product1=join('_',split(" ", $product));
				print transF ">".$tag."|".$product1."\n".$feat_object->spliced_seq->seq."\n"
			}
		}
	}
	close(transF);
}

sub extract_features{
	
	#open CDF,">".$gbk_file."_tRNACDs.txt";
	#open RNAF,">".$gbk_file."_geneDescription.txt";
	open transF,">".$gbk_file."_transcripts.txt";
	
	my %gene=();
	
	while(my $seq_object=$inseq->next_seq){
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			
			if ($feat_object->has_tag("locus_tag")){
				my $tag=($feat_object->get_tag_values("locus_tag"))[0];
				my $gene_name="";
				if($feat_object->has_tag("gene")){
					$gene_name=($feat_object->get_tag_values("gene"))[0];
					$gene{$tag}=$tag."\t".$gene_name."\t".length($feat_object->spliced_seq->seq)."\n";;
				}
				if($gene_name eq ""){
					if(not exists $gene{$tag}){
						$gene{$tag}=$tag."\t".$tag."\t".length($feat_object->spliced_seq->seq)."\n";
					}
				}
			}

			if(($feat_object->primary_tag =~  m/misc_RNA/) || ($feat_object->primary_tag =~  m/rRNA/)){
		        my $tag=($feat_object->get_tag_values("locus_tag"))[0];
		        my $product="N/A";
		        if($feat_object->has_tag("gene")){
		        	$product=($feat_object->get_tag_values("gene"))[0];
		        }elsif($feat_object->has_tag("product")){
		        	$product=($feat_object->get_tag_values("product"))[0];
		        }
				#my $gene_seq=$feat_object->spliced_seq->seq;
				#print RNAF $tag."\t".$product."\n";#|".length($gene_seq)."bp\n".$gene_seq."\n";
				#print transF $tag."\t".length($gene_seq)."\n";
			}
			if($feat_object->primary_tag =~  m/gene/){
				if($feat_object->has_tag("pseudo")){
					my $tag=($feat_object->get_tag_values("locus_tag"))[0];
					#print RNAF $tag."\tpseudo\n";
					print transF $tag."\t".length($feat_object->spliced_seq->seq)."\tpseudo\n";
				}
				if(not ($feat_object->has_tag("gene"))){
					#print transF ($feat_object->get_tag_values("locus_tag"))[0]."\t".length($feat_object->spliced_seq->seq)."\n";
				}
			}
			if(($feat_object->primary_tag =~  m/CDS/)||($feat_object->primary_tag =~  m/RNA/)){
				my $tag=($feat_object->get_tag_values("locus_tag"))[0];
				my $product=($feat_object->get_tag_values("product"))[0];
				#print RNAF $tag."\t".$product."\n";
				print transF $tag."\t".length($feat_object->spliced_seq->seq)."\t".$product."\n";
			}
		}
	}
	#close(RNAF);
	for(keys %gene){
		#print transF $_."\t".$gene{$_}."\n";
	}
	close(transF);
	
	
=pod	
    open RNAF,">".$gbk_file."_rRNA.fasta";
    
	while (my $seq_object = $inseq->next_seq) {
 		for my $feat_object ($seq_object->get_SeqFeatures) { 
 			if($feat_object->has_tag("locus_tag")){
 				print $seq_object->spliced_seq->seq."\n";
 			}

			if($feat_object->primary_tag =~ "rRNA"){
				$gene="N/A";
				if($feat_object->has_tag("gene")){
					$gene=($feat_object->get_tag_values("gene"))[0];
				}
				$tag=($feat_object->get_tag_values("locus_tag"))[0];
				$product=($feat_object->get_tag_values("product"))[0];
				my $gene_seq=$feat_object->spliced_seq->seq;
				print RNAF ">".$tag."\n".$gene_seq."\n";
			}

 		}
	    
	}
	close(RNAF);
=cut
}



sub IMG_gbk_gff_features{
	
    open GFF,">gff_feature.txt";
	while (my $seq_object = $inseq->next_seq) {
 		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if(($feat_object->primary_tag =~ "CDS") || ($feat_object->primary_tag =~ "tRNA")){
				$tag=($feat_object->get_tag_values("locus_tag"))[0];
				$product=($feat_object->get_tag_values("product"))[0];	
				print GFF $tag."\t".$tag."\t".$product."\n";
			}elsif($feat_object->primary_tag =~ "RNA"){
				$gene="N/A";
				if($feat_object->has_tag("gene")){
					$gene=($feat_object->get_tag_values("gene"))[0];
				}
				$tag=($feat_object->get_tag_values("locus_tag"))[0];
				$product=($feat_object->get_tag_values("product"))[0];
				print GFF $gene."\t".$tag."\t".$product."\n";
			}elsif($feat_object->primary_tag =~ "gene"){
				
				if($feat_object->has_tag("pseudo")){
					$product="N/A";
					$tag=($feat_object->get_tag_values("locus_tag"))[0];
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}elsif($feat_object->has_tag("function")){
						$product=($feat_object->get_tag_values("function"))[0];
					}
					print GFF $tag."\t".$tag."\tpseudo:".$product."\n";
				}
			}
 		}
	    
	}
	close(GFF);
}


sub IMG_gbk_to_gff{
	
	my ($gbk_file)	= @_;
    open GFF,">".$gbk_file."_img_parse.gff";
    my %tag_list=();
	while (my $seq_object = $inseq->next_seq) {
		my $genome_id=$seq_object -> id;
		
 		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if(($feat_object->has_tag("locus_tag"))){
				my $start = $feat_object->location->start;
				my $end = $feat_object->location->end;
				my $direc = ($feat_object->strand);
				$tag=($feat_object->get_tag_values("locus_tag"))[0];
				if(exists $tag_list{$tag}){
					next;
				}else{
					$tag_list{$tag}=1;
					my $strand="+";
					if($direc==-1){
						$strand="-";
					}
					my $line=$genome_id."\tRefSeq\tCDS\t".$start."\t".$end."\t.\t".$strand."\t0\ttranscript_id \"".$tag."\"; gene_id \"".$tag."\"; gene_name \"".$tag."\";\n";
					print GFF $line;
				}
			}
 		}
	    
	}
	close(GFF);
	
}


sub genbank_to_fasta{
	
	open FastaFile,">",$gbk_file."_covert.fasta";
	while (my $seq_object = $inseq->next_seq) {
 		my $str_temp=$seq_object->seq;
		$contig_id=$seq_object->id;
		my $total_len=length($str_temp);
	    my $start=0;
	    print FastaFile ">".$contig_id."\n";
	    while(1){
	    	print FastaFile substr($str_temp,$start,60)."\n";
	    	$start=$start+60;
	    	if($start>=$total_len){
	    		last;
	    	}
	    }
	}
	close(FastaFile);
}
sub sub_sequence{
	

	#my $base_dir="/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/ncbi_publication_split_gbkFile";
	# write each entry in the input file to the output file
	my $min_size=500;
	my $seq_out = Bio::SeqIO->new(
    	-file   => ">".$gbk_file."_bp".$min_size.".gbk",
        -format => "genbank",
    );
	
	while (my $seq_object = $inseq->next_seq) {

        if(length($seq_object->seq)>$min_size){
        	$seq_out->write_seq($seq_object);
        }
	    
	}
}

sub split_each_sequence{
	

	#my $base_dir="/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/ncbi_publication_split_gbkFile";
	# write each entry in the input file to the output file
	
	while (my $seq_object = $inseq->next_seq) {
		my $seq_out = Bio::SeqIO->new(
                              -file   => ">"."split_gbkFile/".$seq_object->id.".gbk",
                              -format => "genbank",
                              );
	    $seq_out->write_seq($seq_object);
	}
}

sub split_img_scaffold_annotation{
	
	my $base_dir="/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication";
 	
	open NCBItblFile,">".$base_dir."/gbk_dir/putida_contigs.tbl";
	open NCBIfsaFile,">".$base_dir."/gbk_dir/putida_contigs.fsa";
	open BadProStartFile, "<".$base_dir."/BadProteinStart.id";
	my %BPId=();
	while(<BadProStartFile>){
		chomp($_);
		$BPId{$_}=1;
	}
	close(BadProStartFile);
	open NoStopFile,"<".$base_dir."/NoStopError.id";
	my %NoStopID=();
	while(<NoStopFile>){
		chomp($_);
		$NoStopID{$_}=1;
	}
	#Coordinates of Ns
	my $increase=1;
	my $strSource=" [organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada] [collection-date=Jun-2009] [isolation-source=waste water]";
	my $contig_id="";
	my $count=1;
	my %tbl_line=();
	while(my $seq_object=$inseq->next_seq){
		my $str_temp=$seq_object->seq;
		#print "total length-".length($str_temp)."\n";
		my $scaffold_id=$seq_object->id;
		%ctg_seq=();
		my $start1=0;
		my $total_start=0;
		my $total_length=0;
		if($increase==417){
			#print "pause";
		}
		#open Ctg_CoordFile, ">".$base_dir."/gbk_dir/contig_coord.txt";
		
		while($str_temp  =~ m/(N{100,})/){
			if($increase==11){
				#print "Pause";
			}
			my $find_str=$1;
			#print length($1)."\n";
			$total_length=$total_length+length($1);
			my $result = index($str_temp,$find_str);
			#print "result.".$result."\n";
			my $pre_str=substr($str_temp,0, $result);
			#print "pre_str.".length($pre_str)."\n";
			#print "Ns".length($1)."\n";
			$total_length=$total_length+length($pre_str);
			my $index = sprintf("%05d", $increase);
			$ctg_seq{"contig".$index}="".($start1+1).":".($start1+length($pre_str));
			#print Ctg_CoordFile "contig".$increase."\t".($start1+1)."\t".($start1+length($pre_str))."\t".length($pre_str)."\n";
			$start1=$total_length;
			#print "start.".$start1."\n";
			$str_temp=substr($str_temp,length($pre_str)+length($find_str),);
			print NCBIfsaFile ">contig".$index.$strSource."\n";
			my $index1=0;
			while(($index1+50)<length($pre_str)){
				print NCBIfsaFile substr($pre_str,$index1,50)."\n";
				$index1=$index1+50;
			}
			print NCBIfsaFile substr($pre_str,$index1,)."\n";
			$increase=$increase+1;
			if($increase==12){
				print "Pause";
			}
		}
		my $index = sprintf("%05d", $increase);
		$ctg_seq{"contig".$index}="".($start1+1).":".($start1+length($str_temp));
		print NCBIfsaFile ">contig".$index.$strSource."\n".$str_temp."\n";
		$increase=$increase+1;
		$total_length=$total_length+length($str_temp);
		#print "total_length-".$total_length."\n";
		#close(Ctg_CoordFile);
		#open Gene_CoordFile, ">".$base_dir."/gbk_dir/gene_coordinate.txt";
	
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if($feat_object->primary_tag eq "misc_feature"){
				next;
			}
			if(not $feat_object->has_tag("locus_tag")){
				next;
			}
			if (($feat_object->primary_tag eq "gene") && (not $feat_object->has_tag("pseudo"))){ 
				next;	
			}
			my $start=0;
			my $end=0;
			my $tag="";
			my $img_id="";
			my $product="";
			my $line="";
			my $num_contig="";
			my $start_codon=1;
			my $left=0;
			my $right=0;
			my $direc=0;
			if($feat_object->has_tag("locus_tag")){
				$start = $feat_object->location->start;       
	        	$end = $feat_object->location->end;
	        	$tag=($feat_object->get_tag_values("locus_tag"))[0];
	        	if($tag eq "or1278a"){
	        		#print "pause";
	        	}
	        	my $note=($feat_object->get_tag_values("note"))[0];
	        	$img_id=(split('=',$note))[1];
	        	$direc = ($feat_object->strand);
	        	#print Gene_CoordFile $tag."\t".$start."\t".$end."\t".$direc."\n";
	 			foreach my $temp_id (keys %ctg_seq){
	 				my @arr=split(':',$ctg_seq{$temp_id});
					my $ctg_len=$arr[1]-$arr[0]+1;
					if (($start<=$arr[1]) && ($start >= $arr[0]) && ($end <= $arr[1]) && ($end>= $arr[0])){
						if(($start-$arr[0]+1<4) && ($direc ==1)){
							$start_codon=$start-$arr[0]+1;
							$left=1;					
						}
						if((($arr[1]-$end+1) < 4) && ($direc == 1)){
							$right=1;
						}
						if((($start-$arr[0]+1) < 4) && ($direc == -1)){
							$right=1;
						}
						if((($arr[1]-$end+1) < 4) && ($direc == -1)){
							$left=1;
							$start_codon=$arr[1]-$end+1;
						}	
						$start=$start-$arr[0]+1;
						$end=$end-$arr[0]+1;				
	   				    if($direc==-1){
	 				    	($start,$end)=($end,$start);
	 				    }
	 				    $contig_id=$temp_id;
	 					last;
					}   	
		    	}
		    	$line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				if($right==1){
					$line=$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if($left==1){
					$line="<".$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($right==1)){
					$line="<".$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
			}else{
				next;
			}
			if (($feat_object->primary_tag eq "gene") && ($feat_object->has_tag("pseudo"))){ 
				$line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tpseudo\n";
			}
		    if ($feat_object->primary_tag eq "CDS"){
				if($feat_object->has_tag("product")){
					$product=($feat_object->get_tag_values("product"))[0];
				}else{
					$product="N/A";
				}
				$product =~ s/\.//g;
				#transl_table    11
				my $product_part=$start."\t".$end."\tCDS\n";
				if(($right==1)&&($direc==-1)){
					$product_part=$start."\t>".$end."\tCDS\n";
				}
				if(($left==1)&&($direc==-1)){
					$product_part="<".$start."\t".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($right==1)&&($direc==1)){
					$product_part=$start."\t>".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($direc==1)){
					$product_part="<".$start."\t".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}	
				if(($left==1)&&($$right==1)){
					$product_part="<".$start."\t>".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}					
				#$line=$line.$product_part."\t\t\tcodon_start\t".$start_codon."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|jgi|".$tag."\n\t\t\ttransl_table\t11\n";
		    	if($start_codon==2||$start_codon==3){
		    		$line=$line.$product_part."\t\t\tcodon_start\t".$start_codon."\n\t\t\tnote\tpartial gene\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n\t\t\ttransl_table\t11\n";
		    	}else{
					$line=$line.$product_part."\t\t\tcodon_start\t".$start_codon."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n\t\t\ttransl_table\t11\n";
		    	}
		    }
		    if ($feat_object->primary_tag eq "rRNA"){
				if($feat_object->has_tag("product")){
					$product=($feat_object->get_tag_values("product"))[0];
				}else{
					$product="N/A";
				}
				$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
		    }
		    if ($feat_object->primary_tag eq "tRNA"){
				if($feat_object->has_tag("product")){
					$product=($feat_object->get_tag_values("product"))[0];
				}else{
					$product="N/A";
				}
				$line=$line.$start."\t".$end."\ttRNA\n"."\t\t\tproduct\t".$product."\n";
		    }
		    #note    IMG gene_oid=2510045889
		    
		    if(exists $NoStopID{$tag}){
		    	my $gene_seq=$feat_object->spliced_seq->seq;
				my $stop=0;
				my $index2=0;
				while($index2 <= length($gene_seq)-3){
					my $three=substr($gene_seq,$index2,3);
					if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						print "Scaffold_id ".$scaffold_id." Locus tag ".$tag." product - ".$product. " has stop codon - ".$three." at position ".($index+1)."\n";
						$stop =1;
						last;
					}
					$index2 = $index2+3;
				}
				if(($stop==0) && (($right!=1 && $direc==1) || ($left!=1 && $direc ==-1))){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}
					if($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
					print "scaffold - ".$scaffold_id." locus_tag - ".$tag." product - ".$product." - no stop codon - gene length ".length($gene_seq)." end 10 bases ".substr($gene_seq,length($gene_seq)-11,)."\n";
				}
		    	
		    }
		    
		    if(exists $BPId{$tag}){
		    	my $gene_seq=$feat_object->spliced_seq->seq;
		    	my $first3=substr($gene_seq,0,3);
		    	if ((($left!=1 && $direc==1) || ($right!=1 && $direc ==-1)) && (not $first3 eq "ATG") && (not $first3 eq "GTG") && (not $first3 eq "TTG") && (not $first3 eq "ATT") && 
		    	    (not $first3 eq "CTG") && ($first3 ne "ATA") && ($first3 ne "ATC")){
		    		if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}
					if($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
		    		#print "scaffold_id ".$scaffold_id." locus tag - ".$tag." product - ".$product."\tfirst 3 bases - ".substr($gene_seq,0,3)."\n";
		    	}
		    	
		    }
		    
		    if(not $feat_object->primary_tag =~  m/RNA/){
		    	my $gene_seq=$feat_object->spliced_seq->seq;
				my $first3=substr($gene_seq,0,3);
				if ((($left!=1 && $direc==1) || ($right!=1 && $direc ==-1)) && (not $first3 eq "ATG") && (not $first3 eq "GTG") && (not $first3 eq "TTG") && (not $first3 eq "ATT") && (not $first3 eq "CTG")){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						#print "scaffold - ".$scaffold_id." locus_tag - ".$tag." not start codon - product - ".$product." length ".length($gene_seq)." ".substr($gene_seq,0,10)."\n";

				}	
		    }
		    
		    if(not $feat_object->primary_tag =~  m/RNA/){
		    	my $gene_seq=$feat_object->spliced_seq->seq;
				my $stop=0;
				my $index2=0;
				while($index2 <= length($gene_seq)-3){
					my $three=substr($gene_seq,$index2,3);
					if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
						$stop =1;
						last;
					}
					$index2 = $index2+3;
				}
				if(($stop==0) && (($right!=1 && $direc==1) || ($left!=1 && $direc ==-1))){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}
					if($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
					#print "scaffold - ".$scaffold_id." locus_tag - ".$tag." product - ".$product." - no stop codon - gene length ".length($gene_seq)." end 10 bases ".substr($gene_seq,length($gene_seq)-11,)."\n";
				}
		    }
		    #$line=$line."\n";
		    #print "contig_id - ",$contig_id."\n";
		    $tbl_line{$contig_id}{$count}=$line;
		    $count=$count+1;
		}

		#print "abc";
	}
	#close(Gene_CoordFile);
	
	my @ctg_ids=(keys %tbl_line);
	@ctg_ids=sort(@ctg_ids);
	foreach my $id (@ctg_ids){
		print NCBItblFile ">Feature ".$id."\n";
		foreach my $gid (keys $tbl_line{$id}){
			print NCBItblFile $tbl_line{$id}{$gid};
		}
	}
	close(NCBItblFile);
	close(NCBIfsaFile);
}

sub similar
{
        my($v0,$v1)=@_;
        $v0=~tr/a-z/A-Z/;
        $v1=~tr/a-z/A-Z/;
        
        if (index($v0,".")>-1) { return 0; }
        if (index($v1,".")>-1) { return 0; }
        
        if (index($v0,"HYPO")>-1 || index($v0,"UNCHAR")>-1) { return 0; }
        if (index($v1,"HYPO")>-1 || index($v1,"UNCHAR")>-1) { return 0; }
        
        
        if ($v0 eq $v1) { return 1; }
}


sub overlap_gene{
	
	open TabF, "<ct_contigs.tbl";
	my %gene_h=();
	my $cur_contig="";
	my $cur_start=1000000;
	my $cur_end=1000000;
	my $cur_loc="";
	open OverlapF, ">overlap_gene.txt";
	#open Overlap_2frame_F, ">overlap_gene_same_product_2frame.txt";
	foreach $line (<TabF>){
		chomp($line);
		my $comp=0;
		my $partial=0;
		
		if(not($line =~ ">Feature")){
			if($line =~ ">"){
				$partial=5
			}elsif($line =~ "<"){
				$partial=3
			}
			$line =~ s/<//g;
			$line =~ s/>//g;			
		}
		if($line =~ ">Feature"){
			#print OverlapF $line."\n";
			$cur_contig=$line;
			#$gene_h{$cur_contig}=1;
		}
		if($line =~ "\tgene"){
			my @arr=split("\t",$line);
			$cur_start=$arr[0];
			$cur_end=$arr[1];
		}
		if($line =~ "locus_tag"){
			my @arr=split("\t",$line);
			$gene_h{$cur_contig}{$arr[4]}{"start"}=$cur_start;
			$gene_h{$cur_contig}{$arr[4]}{"end"}=$cur_end;
			
			$cur_loc=$arr[4];
		}
=pod	
		if($line =~ "codon_start"){
			my @arr=split("\t",$line);
			$gene_h{$cur_contig}{$cur_loc}{"codon_start"}=$arr[4];
		}
=cut
		
		my @arr=split("\t",$line);
		if ($arr[3] eq "product"){
			$gene_h{$cur_contig}{$cur_loc}{"product"}=$arr[4];
		}
				
	}
	print "pause";
	my $contig1="";
	my $contig2="";
	my %overlap_h=();
	my $cur_contig="";
	for my $contig (keys %gene_h){
		
		if($contig =~ ">Feature"){
			
            #print $contig."\n";
			for my $gene1 (keys $gene_h{$contig}){
				my $start1=int($gene_h{$contig}{$gene1}{"start"});
				my $end1=int($gene_h{$contig}{$gene1}{"end"});
				#my $condon_start1=$gene_h{$contig}{$gene1}{"codon_start"};
				for my $gene2 (keys $gene_h{$contig}){
					if(($gene1 ne $gene2)){
#					   && (exists $gene_h{$contig}{$gene1}{"codon_start"}) && (exists $gene_h{$contig}{$gene1}{"codon_start"})){
						my $start2=int($gene_h{$contig}{$gene2}{"start"});
						my $end2=int($gene_h{$contig}{$gene2}{"end"});
						my $max1= ($start1,$end1)[$start1<$end1];
						my $min1= ($start1,$end1)[$start1>$end1];
						my $max2= ($start2,$end2)[$start2<$end2];
						my $min2= ($start2,$end2)[$start2>$end2];	
						if((($min2>$min1) && ($min2<$max1)) || (($max2>$min1) && ($max2<$max1))){
							if ((exists $gene_h{$contig}{$gene1}{"product"}) && (exists $gene_h{$contig}{$gene2}{"product"})){
								$product1=$gene_h{$contig}{$gene1}{"product"};
								$product2=$gene_h{$contig}{$gene2}{"product"};
								if(&similar($product1,$product2)==1){
									if( not ((exists $overlap_h{$gene1.$gene2}) || (exists $overlap_h{$gene2.$gene1}))){
										if($cur_contig ne $contig){
											print OverlapF $contig."\n";
											$cur_contig=$contig;
										}
										print OverlapF $gene1.":".$start1." - ".$end1." ".$product1."\t".$gene2.":".$start2." - ".$end2." ".$product2."\n";
										$overlap_h{$gene1.$gene2}=1;
										$overlap_h{$gene2.$gene1}=1;
									}
								}								
							}

							
						}			
=pod
						if(((($start1-$end1)*($start2-$end2)))<0){
							
							next;
						}else{

							if($start1<$end1){
								if((($start2>$start1) && ($start2<$end1)) || (($end2>$start1) && ($end2<$end1))){
									if(abs(($start1-$codon_start1+1)-($start2-$codon_start2+1))%3==0){
										if( not ((exists $overlap_h{$gene1.$gene2}) || (exists $overlap_h{$gene2.$gene1}))){
											print OverlapF $gene1." ".$condon_start1." ".$start1." ".$end1." ".$gene2." ".$condon_start1." ".$start2." ".$end2."\n";
										}
									    $overlap_h{$gene1.$gene2}=1;
									    $overlap_h{$gene2.$gene1}=1;
									}
								}
							}else{
								if((($start2>$end1) && ($start2<$start1)) || (($end2>$end1) && ($end2<$start1))){
									if(abs(($end1-$codon_start1+1)-($end2-$codon_start2+1))%3==0){
										if( not ((exists $overlap_h{$gene1.$gene2}) || (exists $overlap_h{$gene2.$gene1}))){
											print OverlapF $gene1." ".$condon_start1." ".$start1." ".$end1." ".$gene2." ".$condon_start1." ".$start2." ".$end2."\n";
										}
									    $overlap_h{$gene1.$gene2}=1;
									    $overlap_h{$gene2.$gene1}=1;
									}
								}
							}
							
						}
=cut

						
					}
				} 
			}
			
		}
	}
	close(OverlapF);
	
}

sub add_img_dbLink{
	
	my $img_tbl_file="/local/workspace01/zhangju/ct_annotation/gbk_dir/2514885040_CTER.tbl";
    my %protein_db_link=&img_db_link($img_tbl_file);
    open TabF, "<ct_contigs.tbl";
    open NewTabF,">ct_contigs_dbLink.tbl";
    my $old_line="";
    my $cur_tag="";
    foreach $line(<TabF>){
    	chomp($line);
    	my @arr=split("\t",$line);
        if($cur_tag eq "or0025"){
        	print "pause";
        }
    	if($line =~ ">Feature"){
    		if($cur_tag ne ""){
			    if(exists $protein_db_link{$cur_tag} ){
			    	print NewTabF $protein_db_link{$cur_tag};
			    }
    			$cur_tag="";
    		}
     		print NewTabF $line."\n";
     		next;
    	}
    	if($line =~ "locus_tag"){
    		$cur_tag=$arr[4];
    		print NewTabF $line."\n";
    		if($cur_tag eq "or0025"){
        		print "pause";
        	}
    		next;
    	}
    	if(($line =~ "\tgene") && ($cur_tag ne "") ){
		    if(exists $protein_db_link{$cur_tag} ){
		    	print NewTabF $protein_db_link{$cur_tag};
		    }
    		print NewTabF $line."\n";
    		$cur_tag="";
    		next;
    	}
    	if($line ne ""){
    		print NewTabF $line."\n";
    	}
    	
    }
    if(exists $protein_db_link{$cur_tag} ){
    	print NewTabF $protein_db_link{$cur_tag};
    }
}

sub Cintenstinale_ncbi_publication{
	
	my $base_dir="/local/workspace01/zhangju/C.URNW_annotation/ncbi_publication";
	open NCBItblFile,">".$base_dir."/gbk_dir/contigs.tbl";
	open NCBIfsaFile,">".$base_dir."/gbk_dir/contigs.fsa";
	open InternalStopFile, ">",$base_dir."/QC_report/internalStop.txt";
	open StartFile,">",$base_dir."/QC_report/no_startCodon.txt";
	open StopFile, ">",$base_dir."/QC_report/no_stopCodon.txt";
	open NotMultipleFile, ">".$base_dir."/QC_report/not_multipleCodon.txt";

	
	#Coordinates of Ns
	my $increase=1;
	#my $strSource=" [organism=Clostridium termitidis CT 1112, DSM 5398] [strain=CT1112] [country=Canada] [collection-date=1992] [isolation-source=hindgut, termite Nasutitermes lujae]";
	my $strSource=" [organism=Clostridium intestinale strain URNW] [strain=URNW] [country=Canada] [collection-date=2008] [isolation-source=contaminated cellobiose stock solution]";
	my $contig_id="";
	my $count=1;
	my %tbl_line=();
	while(my $seq_object=$inseq->next_seq){
		my $str_temp=$seq_object->seq;
		my $contig_len=length($str_temp);
		$contig_id=$seq_object->id;
		my $start1=0;
		my $total_length=0;
		print NCBIfsaFile ">".$contig_id.$strSource."\n";
		my $index1=0;
		while(($index1+50)<length($str_temp)){
			print NCBIfsaFile substr($str_temp,$index1,50)."\n";
			$index1=$index1+50;
		}
		print NCBIfsaFile substr($str_temp,$index1,)."\n";		
		$total_length=length($str_temp);
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if($feat_object->primary_tag eq "misc_feature"){
				next;
			}
			if(not $feat_object->has_tag("locus_tag")){
				next;
			}
			if (($feat_object->primary_tag eq "gene") && (not $feat_object->has_tag("pseudo"))){ 
				next;	
			}
			my $start=0;
			my $end=0;
			my $tag="";
			my $img_id="";
			my $product="";
			my $line="";
			my $num_contig="";
			my $start_codon=1;
			my $codon_start=1;
			my $left=0;
			my $right=0;
			my $direc=0;
			my @note_values=();
			my $ribo_slippage=0;
			if($feat_object->has_tag("locus_tag")){
				if($feat_object->has_tag("note")){
					@note_values=$feat_object->get_tag_values("note");
					for my $note_value(@note_values){
						if($note_value eq "ribosomal slippage"){
							$ribo_slippage=1;
						}
					}
				}
				if($feat_object->has_tag("ribosomal_slippage")){
					$ribo_slippage=1;
				}

				$start = $feat_object->location->start;       
	        	$end = $feat_object->location->end;
	        	$tag=($feat_object->get_tag_values("locus_tag"))[0];
	        	if($tag eq "or1938"){
	        		print "pause";
	        	}
	        	if($feat_object->has_tag("codon_start")){
	    			$codon_start=($feat_object->get_tag_values("codon_start"))[0];
	    	    }
	        	#my $note=($feat_object->get_tag_values("note"))[0];
	        	#$img_id=(split('=',$note))[1];
	        	$direc = ($feat_object->strand);
	        	if($tag eq "or5351"){
	        		#print "pause";
	        	}
	        	my $edge =0;
				if(($start<4) && ($direc ==1)){
					$start_codon=$start;
					$edge=1;
					$left=1;					
				}
				if((($total_length-$end) < 4) && ($direc == 1)){
					$right=1;
					$edge=1;
				}
				if(($start < 4) && ($direc == -1)){
					$edge=1;
					$right=1;
				}
				if((($total_length-$end) < 4) && ($direc == -1)){
					$edge=1;
					$left=1;
					$start_codon=(($end-$start+1) % 3)+1;
				}	
				if($tag eq "or1938"){
					print "pause";
				}
				#no start codon
				
				if((not $feat_object->primary_tag =~  m/RNA/) && ($left!=1)){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $first3=substr($gene_seq,$start_codon-1,3);
		 			if ((($left!=1)) && (not $first3 eq "ATG") && (not $first3 eq "GTG") && (not $first3 eq "TTG") && (not $first3 eq "ATT") && 
			    	    (not $first3 eq "CTG") && ($first3 ne "ATA") && ($first3 ne "ATC")){
			    			if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							print StartFile "contig - ".$contig_id." locus_tag - ".$tag." no start codon - product - ".$product." length ".length($gene_seq)." ".substr($gene_seq,0,10)."\n";
			    	}
				}
   				if($direc==-1){
 				   	($start,$end)=($end,$start);
 				}
			    $line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
			    
				#not multiple codons
		        if((not $feat_object->primary_tag =~  m/RNA/) && ($right!=1 && $left!=1)){
				#if((not $feat_object->primary_tag =~  m/RNA/) && ($ribo_slippage!=1)){
					my $gene_seq=$feat_object->spliced_seq->seq;
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}elsif($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
					#not multiple codons (x3bps)
					if(($right!=1 && $left!=1)&&(length($gene_seq)%3!=0)){
						print NotMultipleFile "not multiple codons: contig - ".$contig_id." length - ".length($str_temp)."; locus_tag - ".$tag."(".$start."..".$end.") length - ".length($gene_seq)."; product - ".$product."\n";
					}
 				}
 				#no stop codon
 				if((not $feat_object->primary_tag =~  m/RNA/) && ($right!=1)){
	 			#if(not $feat_object->primary_tag =~  m/RNA/){
			    	my $gene_seq=$feat_object->spliced_seq->seq;
					my $stop=0;
					#my $index4=$codon_start-1;
					#my $index2=0;
					my $three=substr($gene_seq,length($gene_seq)-3,3);
					if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
						$stop =1;
					}
					if(($stop==0) && $right!=1 ){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						print StopFile "contig - ".$contig_id." locus_tag - ".$tag." product - ".$product." - no stop codon - gene length ".length($gene_seq)." end 10 bases ".substr($gene_seq,length($gene_seq)-11,)."\n";
					}
			    }
				#internal stop codon
				if((not $feat_object->primary_tag =~  m/RNA/) && ($ribo_slippage!=1)){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $stop=0;
					my $index4=$codon_start-1;
					while($index4 <= length($gene_seq)-6){
						my $three=substr($gene_seq,$index4,3);
						$stop=0;
						if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
							$stop =1;
							if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							if(($product ne "pseudo") && $stop==1){
								print InternalStopFile "internal stop codon: contig - ".$contig_id."\tlocus_tag - ".$tag."\tposition - ".($index4+1)."\tcodon - ".$three."\tproduct - ".$product."\n";
							}
						}
						$index4 = $index4+3;
					}
	
				}
				if($right==1){
					$line=$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if($left==1){
					$line="<".$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($right==1)){
					$line="<".$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
			
			}else{
				next;
			}

			if (($feat_object->primary_tag eq "gene") && ($feat_object->has_tag("pseudo"))){ 
				$line=$line."\t\t\tpseudo\n";
				if($feat_object->has_tag("gene_desc")){
					my $gene_desc=($feat_object->get_tag_values("gene_desc"))[0];
					$line=$line."\t\t\tgene_desc\t".$gene_desc."\n"
				}
				if($feat_object->has_tag("note")){
					my $note=($feat_object->get_tag_values("note"))[0];
					$line=$line."\t\t\tnote\t".$note."\n"
				}				
				
			}else{
				#$line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";#."\t\t\tpseudo\n";
				my $product_part=$start."\t".$end;#."\tCDS\n";
				if(($right==1)&&($direc==-1)){
					$product_part=$start."\t>".$end;#."\tCDS\n";
				}
				if(($left==1)&&($direc==-1)){
					$product_part="<".$start."\t".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($right==1)&&($direc==1)){
					$product_part=$start."\t>".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($direc==1)){
					$product_part="<".$start."\t".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}	
				if(($left==1)&&($right==1)){
					$product_part="<".$start."\t>".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}					

			    if ($feat_object->primary_tag eq "CDS"){
			    	
			    	my @EC_num=();
			    	my $function="";
					if($feat_object->has_tag("EC_number")){
						@EC_num=$feat_object->get_tag_values("EC_number");
					}
			    	
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}
					$product =~ s/\.//g;
					#transl_table    11
=pod
					my $product_part=$start."\t".$end."\tCDS\n";
					if(($right==1)&&($direc==-1)){
						$product_part=$start."\t>".$end."\tCDS\n";
					}
					if(($left==1)&&($direc==-1)){
						$product_part="<".$start."\t".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
					if(($right==1)&&($direc==1)){
						$product_part=$start."\t>".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
					if(($left==1)&&($direc==1)){
						$product_part="<".$start."\t".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}	
					if(($left==1)&&($right==1)){
						$product_part="<".$start."\t>".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
=cut					
					#$line=$line.$product_part."\t\t\tcodon_start\t".$start_codon."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|jgi|".$tag."\n\t\t\ttransl_table\t11\n";
			    	$product_part=$product_part."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
			    	#\t\t\ttransl_table\t11\n";
			    	for my $e_value(@EC_num){
		    			$product_part=$product_part."\t\t\tEC_number\t".$e_value."\n";
		    		}
			    	my @notes=();
			    	my $note_part="";
			    	if($feat_object->has_tag("note")){
						@notes=($feat_object->get_tag_values("note"));
					}
				
					for my $note_value(@notes){
						if((not ($note_value =~ "IMG reference gene")) && (not ($note_value =~ "IMG gene_oid"))){
							$note_part=$note_part."\t\t\tnote\t".$note_value."\n"
						}
						
					}
			    	if($codon_start==2||$codon_start==3){
			    		$line=$line.$product_part."\t\t\tcodon_start\t".$codon_start."\n\t\t\tnote\tpartial gene\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n";#\t\t\ttransl_table\t11\n";
			    	}else{
						$line=$line.$product_part."\t\t\tcodon_start\t".$codon_start."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n";#\t\t\ttransl_table\t11\n";
			    	}
			    	#$line=$line.$product_part."\t\t\tprotein_id\tgnl|jgi|".$tag."\n";
			    	if($ribo_slippage==1){
			    		$line=$line."\t\t\texception\tribosomal slippage\n"
			    	}
			    	$line=$line."\t\t\ttransl_table\t11\n".$note_part;
			    	if($tag eq "or5024"){
			    		#print "start - ".$start."\tend - ".$end."\tright - ".$right."\tleft - ".$left."\n".substr($str_temp,$start-1,$end-$start+1)."\n";
			    	}
			    	if(exists $protein_db_link{$tag}){
		    			$line=$line.$protein_db_link{$tag}."\n";
		    		}
			    	
			    }
			    if ($feat_object->primary_tag eq "rRNA"){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}
					$product_part=$product_part."\trRNA\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tproduct\t".$product."\n";
					#$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
					$line=$line.$product_part;#."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
	
			    }
	
			    if ($feat_object->primary_tag eq "misc_RNA"){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="misc RNA";
					}
					$product_part=$product_part."\tmisc_RNA\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tnote\t".$product."\n";
					#$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
					$line=$line.$product_part;#."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
	
			    }
			    if ($feat_object->primary_tag eq "tRNA"){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}
					$product_part=$product_part."\ttRNA\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tproduct\t".$product."\n";
					#$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
					$line=$line.$product_part;#."\t".$end."\ttRNA\n"."\t\t\tproduct\t".$product."\n";
					
					#$line=$line.$start."\t".$end."\ttRNA\n"."\t\t\tproduct\t".$product."\n";
			    }
			}
		    #note    IMG gene_oid=2510045889

		    $tbl_line{$contig_id}{$count}=$line;
		    $count=$count+1;
		}
	}
	my @ctg_ids=(keys %tbl_line);
	@ctg_ids=sort(@ctg_ids);
	foreach my $id (@ctg_ids){
		print NCBItblFile ">Feature ".$id."\n";
		foreach my $gid (keys $tbl_line{$id}){
			print NCBItblFile $tbl_line{$id}{$gid};
		}
	}
	close(NCBItblFile);
	close(StopFile);
	close(NCBIfsaFile);
	close(InternalStopFile);
	close(StartFile);
	close(NotMultipleFile);
}

sub noGap_ncbi_publication{
	
	#my $base_dir="/local/workspace01/zhangju/ct_annotation";
	
	my $base_dir="/local/workspace01/zhangju/TumsongensisYS13/ncbi_republication";
	open NCBItblFile,">".$base_dir."/gbk_dir/ct_contigs.tbl";
	open NCBIfsaFile,">".$base_dir."/gbk_dir/ct_contigs.fsa";
	open InternalStopFile, ">",$base_dir."/QC_report/internalStop.txt";
	open StartFile,">",$base_dir."/QC_report/no_startCodon.txt";
	open StopFile, ">",$base_dir."/QC_report/no_stopCodon.txt";
	open NotMultipleFile, ">".$base_dir."/QC_report/not_multipleCodon.txt";
	
	my $img_tbl_file="/local/workspace01/zhangju/ct_annotation/2514885040.tbl";
    my %protein_db_link=&img_db_link($img_tbl_file);
	
	my $IMGGeneFastaFile="/local/workspace01/zhangju/ct_annotation/img_gene.fna";
	my %img_gene=();
	our $gene_fasta_in=Bio::SeqIO->new(-file => $IMGGeneFastaFile,
	                               -format => "Fasta"
	                              );
	while (my $seq = $gene_fasta_in->next_seq) {
		$img_gene{(split(' ', $seq->desc))[0]}=$seq->seq;
	}
	#my $img_tbl_file="/local/workspace01/zhangju/wc1_annotation/IMG_generated_ncbi_publication_files/2512047032.tbl";
    my %protein_db_link=&img_db_link($img_tbl_file);
	
	#Coordinates of Ns
	my $increase=1;
	my $strSource=" [organism=Clostridium termitidis CT 1112, DSM 5398] [strain=CT1112] [country=Canada] [collection-date=1992] [isolation-source=hindgut, termite Nasutitermes lujae]";
	my $contig_id="";
	my $count=1;
	my %tbl_line=();
	while(my $seq_object=$inseq->next_seq){
		my $str_temp=$seq_object->seq;
		my $contig_len=length($str_temp);
		$contig_id=$seq_object->id;
		my $start1=0;
		my $total_length=0;
		print NCBIfsaFile ">".$contig_id.$strSource."\n";
		my $index1=0;
		while(($index1+50)<length($str_temp)){
			print NCBIfsaFile substr($str_temp,$index1,50)."\n";
			$index1=$index1+50;
		}
		print NCBIfsaFile substr($str_temp,$index1,)."\n";		
		$total_length=length($str_temp);
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if($feat_object->primary_tag eq "misc_feature"){
				next;
			}
			if(not $feat_object->has_tag("locus_tag")){
				next;
			}
			if (($feat_object->primary_tag eq "gene") && (not $feat_object->has_tag("pseudo"))){ 
				next;	
			}
			my $start=0;
			my $end=0;
			my $tag="";
			my $img_id="";
			my $product="";
			my $line="";
			my $num_contig="";
			my $start_codon=1;
			my $codon_start=1;
			my $left=0;
			my $right=0;
			my $direc=0;
			my @note_values=();
			my $ribo_slippage=0;
			if($feat_object->has_tag("locus_tag")){
				if($feat_object->has_tag("note")){
					@note_values=$feat_object->get_tag_values("note");
					for my $note_value(@note_values){
						if($note_value eq "ribosomal slippage"){
							$ribo_slippage=1;
						}
					}
				}
				if($feat_object->has_tag("ribosomal_slippage")){
					$ribo_slippage=1;
				}

				$start = $feat_object->location->start;       
	        	$end = $feat_object->location->end;
	        	$tag=($feat_object->get_tag_values("locus_tag"))[0];
	        	if($tag eq "or1938"){
	        		print "pause";
	        	}
	        	if($feat_object->has_tag("codon_start")){
	    			$codon_start=($feat_object->get_tag_values("codon_start"))[0];
	    	    }
	        	#my $note=($feat_object->get_tag_values("note"))[0];
	        	#$img_id=(split('=',$note))[1];
	        	$direc = ($feat_object->strand);
	        	if($tag eq "or5351"){
	        		#print "pause";
	        	}
	        	my $edge =0;
				if(($start<4) && ($direc ==1)){
					$start_codon=$start;
					$edge=1;
					$left=1;					
				}
				if((($total_length-$end) < 4) && ($direc == 1)){
					$right=1;
					$edge=1;
				}
				if(($start < 4) && ($direc == -1)){
					$edge=1;
					$right=1;
				}
				if((($total_length-$end) < 4) && ($direc == -1)){
					$edge=1;
					$left=1;
					$start_codon=(($end-$start+1) % 3)+1;
				}	
				if($tag eq "or1938"){
					print "pause";
				}
				#no start codon
				
				if((not $feat_object->primary_tag =~  m/RNA/) && ($left!=1)){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $first3=substr($gene_seq,$start_codon-1,3);
		 			if ((($left!=1)) && (not $first3 eq "ATG") && (not $first3 eq "GTG") && (not $first3 eq "TTG") && (not $first3 eq "ATT") && 
			    	    (not $first3 eq "CTG") && ($first3 ne "ATA") && ($first3 ne "ATC")){
			    			if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							print StartFile "contig - ".$contig_id." locus_tag - ".$tag." no start codon - product - ".$product." length ".length($gene_seq)." ".substr($gene_seq,0,10)."\n";
			    	}
				}
   				if($direc==-1){
 				   	($start,$end)=($end,$start);
 				}
			    $line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
			    
				#not multiple codons
		        if((not $feat_object->primary_tag =~  m/RNA/) && ($right!=1 && $left!=1)){
				#if((not $feat_object->primary_tag =~  m/RNA/) && ($ribo_slippage!=1)){
					my $gene_seq=$feat_object->spliced_seq->seq;
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}elsif($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
					#not multiple codons (x3bps)
					if(($right!=1 && $left!=1)&&(length($gene_seq)%3!=0)){
						print NotMultipleFile "not multiple codons: contig - ".$contig_id." length - ".length($str_temp)."; locus_tag - ".$tag."(".$start."..".$end.") length - ".length($gene_seq)."; product - ".$product."\n";
					}
 				}
 				#no stop codon
 				if((not $feat_object->primary_tag =~  m/RNA/) && ($right!=1)){
	 			#if(not $feat_object->primary_tag =~  m/RNA/){
			    	my $gene_seq=$feat_object->spliced_seq->seq;
					my $stop=0;
					#my $index4=$codon_start-1;
					#my $index2=0;
					my $three=substr($gene_seq,length($gene_seq)-3,3);
					if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
						$stop =1;
					}
					if(($stop==0) && $right!=1 ){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						print StopFile "contig - ".$contig_id." locus_tag - ".$tag." product - ".$product." - no stop codon - gene length ".length($gene_seq)." end 10 bases ".substr($gene_seq,length($gene_seq)-11,)."\n";
					}
			    }
				#internal stop codon
				if((not $feat_object->primary_tag =~  m/RNA/) && ($ribo_slippage!=1)){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $stop=0;
					my $index4=$codon_start-1;
					while($index4 <= length($gene_seq)-6){
						my $three=substr($gene_seq,$index4,3);
						$stop=0;
						if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
							$stop =1;
							if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							if(($product ne "pseudo") && $stop==1){
								print InternalStopFile "internal stop codon: contig - ".$contig_id."\tlocus_tag - ".$tag."\tposition - ".($index4+1)."\tcodon - ".$three."\tproduct - ".$product."\n";
							}
						}
						$index4 = $index4+3;
					}
	
				}
				if($right==1){
					$line=$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if($left==1){
					$line="<".$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($right==1)){
					$line="<".$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
			
			}else{
				next;
			}

			if($tag eq "or5194"){
				print "pause";
			}
			if (($feat_object->primary_tag eq "gene") && ($feat_object->has_tag("pseudo"))){ 
				$line=$line."\t\t\tpseudo\n";
				if($feat_object->has_tag("gene_desc")){
					my $gene_desc=($feat_object->get_tag_values("gene_desc"))[0];
					$line=$line."\t\t\tgene_desc\t".$gene_desc."\n"
				}
				if($feat_object->has_tag("note")){
					my $note=($feat_object->get_tag_values("note"))[0];
					$line=$line."\t\t\tnote\t".$note."\n"
				}				
				
			}else{
				#$line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";#."\t\t\tpseudo\n";
				my $product_part=$start."\t".$end;#."\tCDS\n";
				if(($right==1)&&($direc==-1)){
					$product_part=$start."\t>".$end;#."\tCDS\n";
				}
				if(($left==1)&&($direc==-1)){
					$product_part="<".$start."\t".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($right==1)&&($direc==1)){
					$product_part=$start."\t>".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($direc==1)){
					$product_part="<".$start."\t".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}	
				if(($left==1)&&($right==1)){
					$product_part="<".$start."\t>".$end;#."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
				}					

			    if ($feat_object->primary_tag eq "CDS"){
			    	
			    	my @EC_num=();
			    	my $function="";
					if($feat_object->has_tag("EC_number")){
						@EC_num=$feat_object->get_tag_values("EC_number");
					}
			    	
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}
					$product =~ s/\.//g;
					#transl_table    11
=pod
					my $product_part=$start."\t".$end."\tCDS\n";
					if(($right==1)&&($direc==-1)){
						$product_part=$start."\t>".$end."\tCDS\n";
					}
					if(($left==1)&&($direc==-1)){
						$product_part="<".$start."\t".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
					if(($right==1)&&($direc==1)){
						$product_part=$start."\t>".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
					if(($left==1)&&($direc==1)){
						$product_part="<".$start."\t".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}	
					if(($left==1)&&($right==1)){
						$product_part="<".$start."\t>".$end."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
=cut					
					#$line=$line.$product_part."\t\t\tcodon_start\t".$start_codon."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|jgi|".$tag."\n\t\t\ttransl_table\t11\n";
			    	$product_part=$product_part."\tCDS\n"."\t\t\tlocus_tag\t".$tag."\n";
			    	#\t\t\ttransl_table\t11\n";
			    	for my $e_value(@EC_num){
		    			$product_part=$product_part."\t\t\tEC_number\t".$e_value."\n";
		    		}
			    	my @notes=();
			    	my $note_part="";
			    	if($feat_object->has_tag("note")){
						@notes=($feat_object->get_tag_values("note"));
					}
				
					for my $note_value(@notes){
						if((not ($note_value =~ "IMG reference gene")) && (not ($note_value =~ "IMG gene_oid"))){
							$note_part=$note_part."\t\t\tnote\t".$note_value."\n"
						}
						
					}

			    	if($codon_start==2||$codon_start==3){
			    		$line=$line.$product_part."\t\t\tcodon_start\t".$codon_start."\n\t\t\tnote\tpartial gene\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n";#\t\t\ttransl_table\t11\n";
			    	}else{
						$line=$line.$product_part."\t\t\tcodon_start\t".$codon_start."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n";#\t\t\ttransl_table\t11\n";
			    	}
			    	#$line=$line.$product_part."\t\t\tprotein_id\tgnl|jgi|".$tag."\n";
			    	if($ribo_slippage==1){
			    		$line=$line."\t\t\texception\tribosomal slippage\n"
			    	}
			    	$line=$line."\t\t\ttransl_table\t11\n".$note_part;
			    	if($tag eq "or5024"){
			    		#print "start - ".$start."\tend - ".$end."\tright - ".$right."\tleft - ".$left."\n".substr($str_temp,$start-1,$end-$start+1)."\n";
			    	}
			    	if(exists $protein_db_link{$tag}){
		    			$line=$line.$protein_db_link{$tag}."\n";
		    		}
			    	
			    }
			    if ($feat_object->primary_tag eq "rRNA"){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}
					$product_part=$product_part."\trRNA\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tproduct\t".$product."\n";
					#$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
					$line=$line.$product_part;#."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
	
			    }
			    if($tag eq "CTER_1493"){
			    	print "pause";
			    }
			    if ($feat_object->primary_tag eq "misc_RNA"){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="misc RNA";
					}
					$product_part=$product_part."\tmisc_RNA\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tnote\t".$product."\n";
					#$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
					$line=$line.$product_part;#."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
	
			    }
			    if ($feat_object->primary_tag eq "tRNA"){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}
					$product_part=$product_part."\ttRNA\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tproduct\t".$product."\n";
					#$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
					$line=$line.$product_part;#."\t".$end."\ttRNA\n"."\t\t\tproduct\t".$product."\n";
					
					#$line=$line.$start."\t".$end."\ttRNA\n"."\t\t\tproduct\t".$product."\n";
			    }
			}
		    #note    IMG gene_oid=2510045889

		    $tbl_line{$contig_id}{$count}=$line;
		    $count=$count+1;
		}
	}
	my @ctg_ids=(keys %tbl_line);
	@ctg_ids=sort(@ctg_ids);
	foreach my $id (@ctg_ids){
		print NCBItblFile ">Feature ".$id."\n";
		foreach my $gid (keys $tbl_line{$id}){
			print NCBItblFile $tbl_line{$id}{$gid};
		}
	}
	close(NCBItblFile);
	close(StopFile);
	close(NCBIfsaFile);
	close(InternalStopFile);
	close(StartFile);
	close(NotMultipleFile);
}

sub img_annotation_check{
	
	my $base_dir="/local/workspace01/zhangju/wc1_annotation";
	#Coordinates of Ns
	my $increase=1;
	my $strSource=" [organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]";
	my $contig_id="";
	my $count=1;
	while(my $seq_object=$inseq->next_seq){
		my $str_temp=$seq_object->seq;
		my $contig_len=length($str_temp);
		$contig_id=$seq_object->id;
		my $start1=0;
		my $total_length=0;
		#print NCBIfsaFile ">".$contig_id.$strSource."\n".$str_temp."\n";
		$total_length=length($str_temp);
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if($feat_object->primary_tag eq "misc_feature"){
				next;
			}
			if(not $feat_object->has_tag("locus_tag")){
				next;
			}
			if (($feat_object->primary_tag eq "gene") && (not $feat_object->has_tag("pseudo"))){ 
				next;	
			}
			my $start=0;
			my $end=0;
			my $tag="";
			my $img_id="";
			my $product="";
			my $line="";
			my $num_contig="";
			my $start_codon=1;
			my $left=0;
			my $right=0;
			my $direc=0;
			if($feat_object->has_tag("locus_tag")){
				$start = $feat_object->location->start;       
	        	$end = $feat_object->location->end;
	        	$tag=($feat_object->get_tag_values("locus_tag"))[0];
	        	$direc = ($feat_object->strand);
	        	my $edge =0;
				if(($start<4) && ($direc ==1)){
					$start_codon=$start;
					$edge=1;
					$left=1;					
				}
				if((($total_length-$end) < 4) && ($direc == 1)){
					$right=1;
					$edge=1;
				}
				if(($start < 4) && ($direc == -1)){
					$edge=1;
					$right=1;
				}
				if((($total_length-$end) < 4) && ($direc == -1)){
					$edge=1;
					$left=1;
					$start_codon=(($end-$start+1) % 3)+1;
				}	
				if(not $feat_object->primary_tag =~  m/RNA/){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $first3=substr($gene_seq,0+$start_codon-1,3);
		 			if ((($left!=1)) && (not $first3 eq "ATG") && (not $first3 eq "GTG") && (not $first3 eq "TTG") && (not $first3 eq "ATT") && 
			    	    (not $first3 eq "CTG") && ($first3 ne "ATA") && ($first3 ne "TTG") && ($first3 ne "ATT")){
			    			if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							print "contig - ".$contig_id." length - ".length($str_temp)."; locus_tag - ".$tag."(".$start."..".$end.") length - ".length($gene_seq).". first10 bp ".substr($gene_seq,0,10)." no start codon - product - ".$product."\n";
			    	}
				}
   				if($direc==-1){
 				   	($start,$end)=($end,$start);
 				}
 				
 				if(not $feat_object->primary_tag =~  m/RNA/){
					my $gene_seq=$feat_object->spliced_seq->seq;
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}elsif($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
					#not multiple codons (x3bps)
					if(($start_codon==1)&&(length($gene_seq)%3!=0)){
						#print "contig - ".$contig_id." length - ".length($str_temp)."; locus_tag - ".$tag."(".$start."..".$end.") length - ".length($gene_seq)."; product - ".$product."\n";
					}
 				}

			    if(not $feat_object->primary_tag =~  m/RNA/){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $last_three=substr($gene_seq,length($gene_seq)-3,3);
					if(($right != 1) && (not (($last_three eq "TAA") || ($last_three eq "TAG") || ($last_three eq "TGA")))){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						#print "contig - ".$contig_id."\tlocus_tag - ".$tag."\tlength ".length($gene_seq)." ".substr($gene_seq,length($gene_seq)-10,)."\tproduct - ".$product."\thas no stop codon: TAA/TAG/TGA\n";
					}
				}
	
				if(not $feat_object->primary_tag =~  m/RNA/){
					if($tag eq "TthWC1_0467"){
						print "pause";
					}
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $stop=0;
				    my $index4=$start_codon-1;
					while($index4 <= length($gene_seq)-3){
						my $three=substr($gene_seq,$index4,3);
						if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
							$stop =1;
							if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							#print "contig - ".$contig_id."\tlocus_tag - ".$tag."\tposition - ".($index+1)."\tcodon - ".$three."\tproduct - ".$product."\n";
							last;
						}
						$index4 = $index4+3;
					}
					if(($stop==0) && ($right != 1)){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						#print "contig - ".$contig_id."\tlocus_tag - ".$tag."\tlength ".length($gene_seq)." ".substr($gene_seq,length($gene_seq)-10,)."\tproduct - ".$product."\thas no stop codon: TAA/TAG/TGA\n";
					}
				}
				if($right==1){
					$line=$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if($left==1){
					$line="<".$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
				if(($left==1)&&($right==1)){
					$line="<".$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
				}
			
			}else{
				next;
			}
		    $count=$count+1;
		}
	}
}

sub img_db_link{
	my $img_tbl_file=$_[0];
	my %protein_db_link=();
	
	open tblFile, "<",$img_tbl_file;
	my $tag="";
	my $add_line="";
	my $last_added="";
	foreach my $line (<tblFile>){
		chomp $line;
		
		if(((split(/\t/,$line))[3]) eq "locus_tag"){
			if(($tag ne "") && ($add_line ne "")){
				$protein_db_link{$tag}=$add_line;
			}
			$tag=(split(/\t/,$line))[4];
			
			$add_line="";
		}
=pod
		if((((split(/\t/,$line))[3]) eq "note") && (not ($line =~ "IMG reference gene") && (not ($line =~ "IMG gene_oid")))){
			if($last_added ne ((split(/\t/,$line))[4])){
				$add_line=$add_line."\t\t\tnote\t".((split(/\t/,$line))[4])."\n";
				$last_added=((split(/\t/,$line))[4]);
			}
		}
=cut
		if(((split(/\t/,$line))[3]) eq "db_xref"){
			$add_line=$add_line."\t\t\tdb_xref\t".((split(/\t/,$line))[4])."\n";	
		}
	}
	return %protein_db_link;
}

sub rewrite_gbk{
	my $gbk_in = Bio::SeqIO->new(-file => ">rewrite_gbk.gbk",-format=>"genbank");
	while(my $seq_object=$inseq->next_seq){
		$gbk_in->write_seq($seq_object);	
	}
	

}

sub concat_gbk_custom{
	

	my $gbk_in = Bio::SeqIO->new(-file => ">concate_gbk.gbk",-format=>"genbank");
	my $total_seq="";
	while(my $seq_object=$inseq->next_seq){
		$total_seq=$total_seq.$seq_object->seq."NNNN";
	}
	
	my $total_seq_obj = Bio::Seq->new(-seq=>$total_seq,
									  -display_id => "Geodeb_concatenated");
	my $inseq1=Bio::SeqIO->new(-file => $gbk_file,
                               -format => "genbank"
                             );
	my $start_index=0;
	while(my $seq_object=$inseq1->next_seq){
		for my $feat_object ($seq_object->get_SeqFeatures){
			my $tag_type=$feat_object->primary_tag;
			if(($tag_type eq "gene") || ($tag_type eq "CDS") || ($tag_type =~ "RNA")){
				my %tag_value=();
				my @tag1=();
				my @value1=();
				my $feat = new Bio::SeqFeature::Generic(-start       => $feat_object->start+$start_index,
                                        -end         => $feat_object->end+$start_index,
                                        -strand      => $feat_object->strand,
                                        -primary_tag => $tag_type);
                                        
			    for my $tag ($feat_object->get_all_tags) {             
			        #print "  tag: ", $tag, "\n"; 
			        $tag_value{$tag}=$feat_object->get_tag_values($tag);  
			        if(($tag ne "start") && ($tag ne "end")){
			        	$feat->add_tag_value($tag,($feat_object->get_tag_values($tag))[0]); 
			        }
	        
			    }       
				
				# create the feature with some data, evidence and a note
				
				$total_seq_obj->add_SeqFeature($feat)
			}
		}
		$start_index=$start_index+length($seq_object->seq)+4
	}

	$gbk_in->write_seq($total_seq_obj);	
}

sub concat_gbk_sequence{
	
	#my $io = Bio::SeqIO->new(-file=>'test.gb',-format=>"genbank");#$gbk_file+'_concatenated.gbk');
	my $io = Bio::SeqIO->new(-format => "genbank", -file => ">".$gbk_file."_concatenated.gbk" );       
	my @seqs=();
	while(my $seq_object=$inseq->next_seq){
		push @seqs, $seq_object;
	}
	Bio::SeqUtils->cat(@seqs);
	my $catseq=$seqs[0];
	$io->write_seq($catseq);
	print "test";
}

sub withGap_ncbi_publication{
	
	#my $gname="putida";
	my $gname="wc1";
	my $base_dir="/local/workspace01/zhangju/wc1_annotation";
	#my $base_dir="/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication";
	#open NCBItblFile,">".$base_dir."/gbk_dir/".$gname."_contigs.tbl";
	open NCBIfsaFile,">".$base_dir."/gbk_dir/".$gname."_contigs.fsa";
	open NCBIGeneFile,">".$base_dir."/gbk_dir/".$gname."_contigs.gene";
	#open StopFile, ">",$base_dir."/".$gname."_internalStop.txt";
	
	open InternalStopFile, ">",$base_dir."/QC_report/internalStop.txt";
	open StartFile,">",$base_dir."/QC_report/no_startCodon.txt";
	open StopFile, ">",$base_dir."/QC_report/no_stopCodon.txt";
	open NotMultipleFile, ">".$base_dir."/QC_report/not_multipleCodon.txt";
	open AGPFile, ">".$base_dir."/gbk_dir/".$gname."_scaffolds.agp";
	my $IMGGeneFastaFile=$base_dir."/img_gene.fna";
	my %img_gene=();
	our $gene_fasta_in=Bio::SeqIO->new(-file => $IMGGeneFastaFile,
	                               -format => "Fasta"
	                              );
	while (my $seq = $gene_fasta_in->next_seq) {
		$img_gene{(split(' ', $seq->desc))[0]}=$seq->seq;
	}
	my $img_tbl_file="/local/workspace01/zhangju/wc1_annotation/IMG_generated_ncbi_publication_files/2512047032.tbl";
    my %protein_db_link=&img_db_link($img_tbl_file);
	#Coordinates of Ns
	my $increase=1;
	#" [organism=Clostridium termitidis CT 1112, DSM 5398] [strain=CT1112] [country=Africa: Congo] [collection-date=1992] [isolation-source=hindgut, termite Nasutitermes lujae]"
	my $strSource=" [organism=Thermoanaerobacter thermohydrosulfuricus] [strain=WC1] [country=Canada] [collection-date=Jul-2006] [isolation-source=woodchip compost]";
	my $contig_id="";
	my $count=1;
	my %tbl_line=();
	my %contig_dna=();
	my %scaffold_map=();
	my $scaffold_index=1;
	while(my $seq_object=$inseq->next_seq){
		my $str_temp=$seq_object->seq;
		if(length($str_temp)<200){
			next;
		}
		
		my $NCBI_scaffoldId = "scaffold".sprintf("%05d", $scaffold_index);
		my $scaffold_id=$seq_object->id;
		$scaffold_map{$NCBI_scaffoldId}=$scaffold_id;
		#print "scaffold mapping - "$NCBI_scaffoldId."\t".$scaffold_id."\n";
		#print "total length-".length($str_temp)."\n";
		if($scaffold_id eq "twc1a_scaffold00042.42"){
			print "pause";
		}
=for comment
		if($scaffold_id eq "PLS46_scaffold00001.1"){
			open SC1File , ">".$base_dir."/scaffold1.fna";
			print SC1File ">".$scaffold_id."\n";
			print SC1File $str_temp."\n";
			close(SC1File);
		}
=cut
		my $scaffold_length=length($str_temp);
		%ctg_seq=();
		my $start1=0;
		my $total_start=0;
		my $total_length=0;

		#open Ctg_CoordFile, ">".$base_dir."/gbk_dir/contig_coord.txt";
=for comment		
		while($str_temp  =~ m/(N{1,})/){
			my $find_str=$1;
			#print length($1)."\n";
			#
			my $result = index($str_temp,$find_str);
			my $pre_str=substr($str_temp,0, $result);
			$str_temp=substr($str_temp,length($pre_str)+length($find_str),);
			#print "scaffold_id ".$scaffold_id." number of Ns ".length($1)." start from position ".($total_length+$result)."\n";
			$total_length=$total_length+length($1)+length($pre_str);
		}
=cut
		
		$str_temp=$seq_object->seq;
		$total_length=0;
		my $ncbi_scaffold_contig_count=1;
		while($str_temp  =~ m/(N{5,})/){
			my $find_str=$1;
			#print length($1)."\n";
			$total_length=$total_length+length($1);
			my $result = index($str_temp,$find_str);
			#print "result.".$result."\n";
			my $pre_str=substr($str_temp,0, $result);
			#print "pre_str.".length($pre_str)."\n";
			#print "Ns".length($1)."\n";
			if(length($pre_str)<200){
				$total_length=$total_length+length($pre_str);
				$start1=total_length;
				print "short than 200 - ".$scaffold_id." contig ".sprintf("%05d", $increase);
				next;
			}
			$total_length=$total_length+length($pre_str);
			my $index = sprintf("%05d", $increase);
			#$ctg_seq{"contig".$index."_".$scaffold_id}="".($start1+1).":".($start1+length($pre_str));
			$ctg_seq{"contig".$index}="".($start1+1).":".($start1+length($pre_str));
			#print Ctg_CoordFile "contig".$increase."\t".($start1+1)."\t".($start1+length($pre_str))."\t".length($pre_str)."\n";
			my $agp_sequence=$NCBI_scaffoldId."\t".($start1+1)."\t".
			    ($start1+length($pre_str))."\t".$ncbi_scaffold_contig_count."\tW\tcontig".$index."\t1\t".(length($pre_str))."\t+\n";
			$agp_sequence=$agp_sequence.$NCBI_scaffoldId."\t".($start1+length($pre_str)+1)."\t".($start1+length($pre_str)+length($find_str))."\t".($ncbi_scaffold_contig_count+1).
			    "\tN\t".length($find_str)."\tfragment\tyes\n";
			print AGPFile $agp_sequence;
			$start1=$total_length;
			#print "start.".$start1."\n";
			$str_temp=substr($str_temp,length($pre_str)+length($find_str),);
			#print NCBIfsaFile ">contig".$index.$strSource."\n";
			#print NCBIfsaFile ">contig".$index."_".$scaffold_id." ".length($pre_str)."bps first10bp ".substr($pre_str,0,10).
		    #  " last10bp ".substr($pre_str,length($pre_str)-10,)." scaffold_id ".$scaffold_id." length ".$scaffold_length."\n";#$strSource."\n".$str_temp."\n";
			print NCBIfsaFile ">contig".$index.$strSource."\n";
			#print "contig mapping - contig".$index." ".$scaffold_id;
			$contig_dna{"contig".$index}=$pre_str;
			my $index1=0;
			while(($index1+50)<length($pre_str)){
				print NCBIfsaFile substr($pre_str,$index1,50)."\n";
				$index1=$index1+50;
			}
			print NCBIfsaFile substr($pre_str,$index1,)."\n";
			$increase=$increase+1;
			$ncbi_scaffold_contig_count=$ncbi_scaffold_contig_count+2;
		}
		my $index = sprintf("%05d", $increase);
		#$ctg_seq{"contig".$index."_".$scaffold_id}="".($start1+1).":".($start1+length($str_temp));
		$ctg_seq{"contig".$index}="".($start1+1).":".($start1+length($str_temp));
		my $agp_sequence=$NCBI_scaffoldId."\t".($start1+1)."\t".
			    ($start1+length($str_temp))."\t".$ncbi_scaffold_contig_count."\tW\tcontig".$index."\t1\t".(length($str_temp))."\t+\n";
		print AGPFile $agp_sequence;

		#print NCBIfsaFile ">contig".$index."_".$scaffold_id." ".length($str_temp)."bps first10bp ".substr($str_temp,0,10).
		#     " last10bp ".substr($str_temp,length($str_temp)-10,)." scaffold_id ".$scaffold_id." length ".$scaffold_length."\n";#$strSource."\n".$str_temp."\n";
		#print NCBIfsaFile ">".$scaffold_id."_contig".$index.$strSource."\n";
		print NCBIfsaFile ">contig".$index.$strSource."\n";
		$contig_dna{"contig".$index}=$str_temp;
		my $index5=0;
		while(($index5+50)<length($str_temp)){
			print NCBIfsaFile substr($str_temp,$index5,50)."\n";
			$index5=$index5+50;
		}
		print NCBIfsaFile substr($str_temp,$index5,)."\n";
		
		$increase=$increase+1;
		$total_length=$total_length+length($str_temp);
		#print "total_length-".$total_length."\n";
		#close(Ctg_CoordFile);
		#open Gene_CoordFile, ">".$base_dir."/gbk_dir/gene_coordinate.txt";
		my %gap_gene=();
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			if($feat_object->primary_tag eq "misc_feature"){
				next;
			}
			if(not $feat_object->has_tag("locus_tag")){
				next;
			}
			if (($feat_object->primary_tag eq "gene") && (not $feat_object->has_tag("pseudo"))){ 
				next;	
			}
			my $start=0;
			my $end=0;
			my $tag="";
			my $img_id="";
			my $product="";
			my $line="";
			my $num_contig="";
			my $start_codon=1;
			my $codon_start=1;
			my $left=0;
			my $right=0;
			my $direc=0;
			if($feat_object->has_tag("locus_tag")){
				$start = $feat_object->location->start;       
	        	$end = $feat_object->location->end;
	        	$tag=($feat_object->get_tag_values("locus_tag"))[0];
	        	if($tag eq "TthWC1_0253"){
	        		print "pause";
	        	}
=pod
	        	my $gene_seq=$feat_object->spliced_seq->seq;
	        	if($gene_seq  =~ m/(N{1,})/){
					my $find_str=$1;
					$gap_gene{$tag}=length($1);
					my $result = index($gene_seq,$find_str);
					#print "result.".$result."\n";
					my $pre_str=substr($gene_seq,0, $result);
					my $pos_str=substr($gene_seq,length($pre_str.$1)-1,);
					print $tag." length ".length($gene_seq)." has Ns of size ".
					  length($1)." left part - ". length($pre_str)." right part - ".length($pos_str)."\n";
				}
=cut
	        	#my $note=($feat_object->get_tag_values("note"))[0];
	        	#$img_id=(split('=',$note))[1];
	        	$direc = ($feat_object->strand);
	        	#print Gene_CoordFile $tag."\t".$start."\t".$end."\t".$direc."\n";
	        	my $in_contig=0;
	 			foreach my $temp_id (keys %ctg_seq){
	 				my @arr=split(':',$ctg_seq{$temp_id});
					my $ctg_len=$arr[1]-$arr[0]+1;
					if (($start<=$arr[1]) && ($start >= $arr[0]) && ($end <= $arr[1]) && ($end>= $arr[0])){
						if(($start-$arr[0]+1<4) && ($direc ==1)){
							$start_codon=$start-$arr[0]+1;
							$left=1;					
						}
						if((($arr[1]-$end+1) < 4) && ($direc == 1)){
							$right=1;
						}
						if((($start-$arr[0]+1) < 4) && ($direc == -1)){
							$right=1;
						}
						if((($arr[1]-$end+1) < 4) && ($direc == -1)){
							$left=1;
							$start_codon=$arr[1]-$end+1;
						}	
						$start=$start-$arr[0]+1;
						$end=$end-$arr[0]+1;	
						#my $split_gene_seq="";	
=for comment
						if($direc==1){
							$split_gene_seq=substr($contig_dna{$temp_id},$start-1,$end-$start+1);
						}else{
							$split_gene_seq=revcom(substr($contig_dna{$temp_id},$start-1,$end-$start+1))->seq;
							
						}
						if($split_gene_seq ne $img_gene{$tag}){
							print $scaffold_id." ".$contig_id." ".$tag." sequence not idential\n"."img_gene ".$img_gene{$tag}
							     ."\nsplit_gene ".$split_gene_seq."\n";
						}
=cut
	   				    if($direc==-1){
	 				    	($start,$end)=($end,$start);
	 				    }
	 				    $contig_id=$temp_id;
						$in_contig=1;
	 					last;
					}   	
		    	}
		    	if($in_contig==1){
			    	$line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
					if($right==1){
						$line=$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
					if($left==1){
						$line="<".$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
					}
					if(($left==1)&&($right==1)){
						$line="<".$start."\t>".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n";
					}		    		
		    	}else{
		    		print "not in contig ".$tag."\n";
		    	}

			}else{
				next;
			}
			#if (($feat_object->primary_tag eq "gene") && ($feat_object->has_tag("pseudo"))){ 
			#	$line=$start."\t".$end."\tgene\n"."\t\t\tlocus_tag\t".$tag."\n"."\t\t\tpseudo\n";
			#}
			
			if (($feat_object->primary_tag eq "gene") && ($feat_object->has_tag("pseudo"))){ 
				$line=$line."\t\t\tpseudo\n";
			}else{
		    #if ($feat_object->primary_tag eq "CDS"){
		    	#/function="glutamate racemase (EC 5.1.1.3)"
		    	if ($feat_object->primary_tag eq "CDS"){
			    	if($feat_object->has_tag("codon_start")){
			    		$codon_start=($feat_object->get_tag_values("codon_start"))[0];
			    	}
			    	my @EC_num=();
			    	my $function="";
					if($feat_object->has_tag("EC_number")){
						@EC_num=$feat_object->get_tag_values("EC_number");
					}
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}else{
						$product="N/A";
					}

					
					$product =~ s/\.//g;
					#transl_table    11
					my $product_part=$start."\t".$end."\tCDS\n";
					if(($left==1)&&($right==1)){
						$product_part="<".$start."\t>".$end."\tCDS\n";#."\t\t\tlocus_tag\t".$tag."\n";
					}else{
						if($right==1){
							$product_part=$start."\t>".$end."\tCDS\n";
						}
						if($left==1){
							$product_part="<".$start."\t".$end."\tCDS\n"
						}
						
					}			
					#$line=$line.$product_part."\t\t\tcodon_start\t".$start_codon."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|jgi|".$tag."\n\t\t\ttransl_table\t11\n";
			    	for my $e_value(@EC_num){
			    		$product_part=$product_part."\t\t\tEC_number\t".$e_value."\n";
			    	}
			    	if($start_codon==2||$start_codon==3){
			    		$line=$line.$product_part."\t\t\tcodon_start\t".$codon_start."\n\t\t\tnote\tpartial gene\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n";#\t\t\ttransl_table\t11\n";
			    	}else{
						$line=$line.$product_part."\t\t\tcodon_start\t".$codon_start."\n\t\t\tproduct\t".$product."\n\t\t\tprotein_id\tgnl|mgcb2|".$tag."\n";#\t\t\ttransl_table\t11\n";
			    	}
			    	if($feat_object->has_tag("inference")){
			    		$line=$line."\t\t\tinference\t".(($feat_object->get_tag_values("product"))[0])."\n";
			    	}
			    	if($feat_object->has_tag("transl_except")){
			    		$line=$line."\t\t\ttransl_except\n";
			    	}
			    	if($feat_object->has_tag("ribosomal_slippage")){
			    		$line=$line."\t\t\texception\tribosomal slippage\n";
			    	}

			    	$line=$line."\t\t\ttransl_table\t11\n";
			    	if($tag eq "TthWC1_0253"){
			    		print "pause";
			    	}
			    	if($feat_object->has_tag("note")){
			    		$str_note=($feat_object->get_tag_values("note"))[0];
			    		@arr_note=split(";",$str_note);
			    		for(@arr_note){
			    			#$str1 = ((split("=",$_))[1]);
			    			my $s_note=$_;
			    			$s_note =~ s/^\s+//; 
			    			$line=$line."\t\t\tnote\t".$s_note."\n";
			    		}
			    	}
	        	   	if(exists $protein_db_link{$tag}){
			    		$line=$line.$protein_db_link{$tag}."\n";
			    	}
		    	}
		    }
		    if ($feat_object->primary_tag eq "rRNA"){
				if($feat_object->has_tag("product")){
					$product=($feat_object->get_tag_values("product"))[0];
				}else{
					$product="N/A";
				}
				$line=$line.$start."\t".$end."\trRNA\n"."\t\t\tproduct\t".$product."\n";
		    }
		    if ($feat_object->primary_tag eq "tRNA"){
				if($feat_object->has_tag("product")){
					$product=($feat_object->get_tag_values("product"))[0];
				}else{
					$product="N/A";
				}
				$line=$line.$start."\t".$end."\ttRNA\n"."\t\t\tproduct\t".$product."\n";
		    }
		    #note    IMG gene_oid=2510045889
		    #no start codon
		    if((not $feat_object->primary_tag =~  m/RNA/) && ($left!=1)){
		    	my $gene_seq=$feat_object->spliced_seq->seq;
				my $first3=substr($gene_seq,0,3);
				if ((not $first3 eq "ATG") && (not $first3 eq "ATA") && (not $first3 eq "ATC") && (not $first3 eq "GTG")
				    && (not $first3 eq "TTG") && (not $first3 eq "ATT") && (not $first3 eq "CTG")){
						if($feat_object->has_tag("product")){
							$product=($feat_object->get_tag_values("product"))[0];
						}
						if($feat_object->has_tag("pseudo")){
							$product="pseudo";
						}
						if($tag eq "or0015"){
							print "pause";
						}
						if($feat_object->has_tag("inference")){
							
						}else{
							print StartFile "no start codon: scaffold - ".$scaffold_id." locus_tag - ".$tag." - product - ".$product." length ".length($gene_seq)." ".substr($gene_seq,0,10)."\n";
						}
		    	}
		    }

		    #not multiple codons
		    if((not $feat_object->primary_tag =~  m/RNA/) && ($right!=1 && $left!=1)){
				my $gene_seq=$feat_object->spliced_seq->seq;
				if($feat_object->has_tag("product")){
					$product=($feat_object->get_tag_values("product"))[0];
				}elsif($feat_object->has_tag("pseudo")){
					$product="pseudo";
				}
				#not multiple codons (x3bps)
				if((length($gene_seq)%3!=0)){
					print notMutipleFile "not multiple codons: ".$scaffold_id." contig - ".$contig_id." length - ".length($str_temp)."; locus_tag - ".$tag."(".$start."..".$end.") length - ".length($gene_seq)."; product - ".$product."\n";
				}
 			}
 			#internal stop codon
 			if(not $feat_object->primary_tag =~  m/RNA/){
					my $gene_seq=$feat_object->spliced_seq->seq;
					my $stop=0;
				    my $index4=$codon_start-1;
					while($index4 <= length($gene_seq)-6){
						my $three=substr($gene_seq,$index4,3);
						$stop=0;
						if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
							$stop =1;
							if($feat_object->has_tag("product")){
								$product=($feat_object->get_tag_values("product"))[0];
							}
				
							if($feat_object->has_tag("pseudo")){
								$product="pseudo";
							}
							if(($product ne "pseudo") && ($stop==1) && (not $feat_object->has_tag("transl_except"))){
								print InternalStopFile "internal stop codon: contig - ".$contig_id."\tlocus_tag - ".$tag."\tcodon_start\t".$codon_start."\tposition - ".($index4+1)."\tcodon - ".$three."\tproduct - ".$product."\n";
								#print InternalStopFile "internal stop codon: contig - ".$contig_id."\tlocus_tag - ".$tag."\tposition - ".($index4+1)."\tcodon - ".$three."\tproduct - ".$product."\n";
							}
						}
						$index4 = $index4+3;
					}
			}

			#stop codon
			if((not $feat_object->primary_tag =~  m/RNA/) && ($right!=1)){
		    	my $gene_seq=$feat_object->spliced_seq->seq;
				my $stop=0;
				my $index2=0;
				my $three=substr($gene_seq,length($gene_seq)-3,3);
				if(($three eq "TAA") || ($three eq "TAG") || ($three eq "TGA")){
					$stop =1;
				}
				if($stop==0){
					if($feat_object->has_tag("product")){
						$product=($feat_object->get_tag_values("product"))[0];
					}
					if($feat_object->has_tag("pseudo")){
						$product="pseudo";
					}
					print StopFile "contig - ".$contig_id." locus_tag - ".$tag." product - ".$product." - no stop codon - gene length ".length($gene_seq)." end 10 bases ".substr($gene_seq,length($gene_seq)-11,)."\n";
				}
			}
		
		    #$line=$line."\n";
		    #print "contig_id - ",$contig_id."\n";
		    #$tbl_line{$scaffold_id."_".$contig_id}{$count}=$line;
		    $tbl_line{$contig_id}{$count}=$line;
		    $count=$count+1;
		}
		$scaffold_index=$scaffold_index+1;
		#print "abc";
	}
	#close(Gene_CoordFile);
	
	my @ctg_ids=(keys %tbl_line);
	@ctg_ids=sort(@ctg_ids);
	foreach my $id (@ctg_ids){
		print NCBItblFile ">Feature ".$id."\n";
		foreach my $gid (keys $tbl_line{$id}){
			my @arr1=split('\n',$tbl_line{$id}{$gid});
			my $start = (split('\t',$arr1[0]))[0];
			my $end = (split('\t',$arr1[0]))[1];
			my $tag=(split('\t',$arr1[1]))[4];
			print NCBIGeneFile $id."\t".$tag."\t".$start."\t".$end."\n";
			#print NCBItblFile $tbl_line{$id}{$gid};
		}
	}
	close(NCBIGeneFile);
	#close(NCBItblFile);
	close(NCBIfsaFile);
	close(StopFile);
	close(AGPFile);
}


sub checking_ncbi_tbl_on_fsa{
	
	my $base_dir="/local/workspace01/zhangju/wc1_annotation/gbk_dir/wc1_contigs";
	my $fsa_file = $base_dir.".fsa";
	my $gene_file=$base_dir.".gene";
	#my $format="Fasta";
	our $inseq=Bio::SeqIO->new(-file => $fsa_file,
	                               -format => "Fasta"
	                              );
	my %fsa_seq=();
	while (my $seq = $inseq->next_seq) {
		$fsa_seq{$seq->id}=$seq->seq;
	}
	my %gene_seq=();
	
	open GF, "<", $gene_file;
	
	while(<GF>){
		chomp($_);
		my $line=$_;
		$line =~ s/<//;
		$line =~ s/>//;
		my @arr=split("\t",$line);
		my $id= $arr[0];
		my $tag=$arr[1];
		my $start=$arr[2];
		my $end=$arr[3];
		my $geneSeq="";
		if($start < $end){
			my $len=$end-$start+1;
			$geneSeq=substr($fsa_seq{$id},$start-1,$len)
		}else{
			my $len=$start-$end+1;
			$f_seq=substr($fsa_seq{$id},$end-1,$len);
			$geneSeq=revcom($f_seq)->seq;
		}
		$gene_seq{$tag}=$geneSeq;
		if($id eq "contig00002_twc1a_scaffold00002.2"){
			print ">".$tag." length ".length($gene_seq{$tag})." first 10bp ".substr($gene_seq{$tag},0,10)." last 10bp ".substr($gene_seq{$tag},length($gene_seq{$tag})-10,10)."\n";
		}
	}
	close(GF);
}

#generating fasta protein of CDs and pseudo genes
sub each_gbk_protein_fasta{
	#my $gbk_file="ct_contigs.gbk";
	my $fasta=$gbk_file."_protein.faa";
	my $fasta_gene=$gbk_file."_gene.faa";
	my $format="genbank";
	#my $seqio_object=Bio::SeqIO->new(-file=>$gbk_file, -format=>$format);
	open(ParseFile,">".$fasta);
	open(ParseFileGene,">".$fasta_gene);
	open(DiscrepF,">protein_seq_discrepancy.txt");
	my $num_cd=0;
	my $num_pseudo=0;
	while(my $seq_object = $inseq->next_seq){
		print DiscrepF "====".$seq_object->display_id()."====\n";
		for my $feat_object ($seq_object->get_SeqFeatures) { 
			#if ($feat_object->primary_tag eq "LOCUS"){
			#($feat_object->has_tag("pseudo")) || 
		    if ($feat_object->primary_tag eq "CDS"){
		    	my $protein_function="";
		    	if($feat_object->has_tag("product")){
		    		$protein_function="|product - ".($feat_object->get_tag_values('product'))[0]
		    	}
		    	if($feat_object->has_tag("note")){
		    		$protein_function=$protein_function."|note - ".($feat_object->get_tag_values('note'))[0]
		    	}		    	
		    	$num_cd=$num_cd+1;
		        my $tag_value=($feat_object->get_tag_values('locus_tag'))[0];
		        my $translate="n/a";
		        if ($feat_object->primary_tag eq "CDS"){
		        	if (not ($feat_object->has_tag("translation"))){
		        		print $tag_value." has no translation\n"
		        	}else{
		        		$translate=($feat_object->get_tag_values('translation'))[0];
		        	}
		        }
	        	my $seq_obj = Bio::Seq->new(-seq => $feat_object->spliced_seq->seq, -alphabet => 'dna' );
	        	my $dna_seq = $feat_object->spliced_seq->seq;
	        	my $annotated_seq="";
	        	if($tag_value eq "CTER_0018"){
	        		print "pause";
	        	}
	        	my $pro = translate($dna_seq,-codontable_id => 11);
	        	my $pro_seq="";
	        	if($pro->seq =~ /\*$/){
	        		$pro_seq=substr($pro->seq,0,length($pro->seq)-1);
	        	}else{
	        		$pro_seq=$pro->seq;
	        	}
	        	if (not ($pro_seq =~ /\*/)){
	        		#$annotated_seq=$pro->seq;
	        		$annotated_seq=$pro_seq;
	        	}else{
	        		$pro= translate(substr($dna_seq,1,),-codontable_id => 11);
	        		if($pro->seq =~ /\*$/){
	        			$pro_seq=substr($pro->seq,0,length($pro->seq)-1);
		        	}else{
		        		$pro_seq=$pro->seq;
		        	}
	        		if (not ($pro_seq =~ /\*/)){
	        			$annotated_seq=$pro_seq;
	        		}else{
	        			$pro= translate(substr($dna_seq,2,),-codontable_id => 11);
	        			if($pro->seq =~ /\*$/){
		        			$pro_seq=substr($pro->seq,0,length($pro->seq)-1);
			        	}else{
			        		$pro_seq=$pro->seq;
			        	}
	        			if (not ($pro_seq =~ /\*/)){
	        				$annotated_seq=$pro_seq;
	        			}
	        		}
	        	}
	        	if($translate ne "n/a"){
	        		print ParseFile ">".$tag_value."\n";#.$protein_function."\n";
	        		print ParseFile $translate."\n";#  
	        		print ParseFileGene ">".$tag_value.$protein_function."\n";
	        		print ParseFileGene $dna_seq."\n";#       
	        	}else{
	        		print $tag_value." has not translation\n";
	        	}
        		if(substr($annotated_seq,1,) ne substr($translate,1,)){
        			print DiscrepF $tag_value." not equal\n";
        			print DiscrepF "translation    - ".$translate."\n";
        			print DiscrepF "dna translated - ".$annotated_seq."\n";
        		}
		    }
		    if ($feat_object->has_tag("pseudo")){
		    	$num_pseudo=$num_pseudo+1;
		        my $tag_value=($feat_object->get_tag_values('locus_tag'))[0];
		
	        	my $dna_seq = $feat_object->spliced_seq->seq;
	        	my $annotated_seq="";
	    
	        	my $pro = translate($dna_seq,-codontable_id => 11);
	        	my $pro_seq="";
	        	if($pro->seq =~ /\*$/){
	        		$pro_seq=substr($pro->seq,0,length($pro->seq)-1);
	        	}else{
	        		$pro_seq=$pro->seq;
	        	}
	        	if (not ($pro_seq =~ /\*/)){
	        		#$annotated_seq=$pro->seq;
	        		$annotated_seq=$pro_seq;
	        	}else{
	        		$pro= translate(substr($dna_seq,1,),-codontable_id => 11);
	        		if($pro->seq =~ /\*$/){
	        			$pro_seq=substr($pro->seq,0,length($pro->seq)-1);
		        	}else{
		        		$pro_seq=$pro->seq;
		        	}
	        		if (not ($pro_seq =~ /\*/)){
	        			$annotated_seq=$pro_seq;
	        		}else{
	        			$pro= translate(substr($dna_seq,2,),-codontable_id => 11);
	        			if($pro->seq =~ /\*$/){
		        			$pro_seq=substr($pro->seq,0,length($pro->seq)-1);
			        	}else{
			        		$pro_seq=$pro->seq;
			        	}
	        			if (not ($pro_seq =~ /\*/)){
	        				$annotated_seq=$pro_seq;
	        			}
	        		}
	        	}
	        	if($annotated_seq ne ""){
	        		print ParseFile ">".$tag_value."\n";
	        		print ParseFile $annotated_seq."\n";#  
	        		print ParseFileGene ">".$tag_value."\n";
	        		print ParseFileGene $dna_seq."\n";#     
	        	}	        	
		    }
		}
	}
	print "number of cds ".$num_cd."\n";
	print "number of pseudos ".$num_pseudo."\n";
	close(ParseFile);
	close(DiscrepF);
	close(ParseFileGene);
}

#&each_gbk_protein_fasta();
#&overlap_gene();
#&add_img_dbLink();
#exit 0;

our $gbk_file = shift;
my $format="genbank";
our $inseq=Bio::SeqIO->new(-file => $gbk_file,
                               -format => $format
                             );
&extract_genes_canola();
#&extract_genes();
#&IMG_gbk_to_gff($gbk_file);                             
#&Cintenstinale_ncbi_publication();                             
#&each_gbk_protein_fasta();
#&concat_gbk_sequence();
#&rewrite_gbk();
#&concat_gbk_custom();
#&extract_CSterc();
#&extract_gd1_genePRIMP();
#&extract_features();
#&IMG_gbk_gff_features();
#&splitScaffold_wc1_annotation();
#&checking_ncbi_tbl_on_fsa();
#&ct_ncbi_publication();
#&withGap_ncbi_publication();
#&noGap_ncbi_publication();
#&split_each_sequence();
#&sub_sequence();
#&overlap_gene();
&genbank_to_fasta();
exit 0;
#&img_annotation_check();
#&ct_ncbi_publication();
#&split_img_scaffold_annotation();