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

use Bio::SeqIO;
use Bio::Seq;

sub open_fasta{
	(my $fasta_file)=@_;
	my $format="Fasta";
	my $inseq=Bio::SeqIO->new(-file => $fasta_file,
                               -format => $format
                              );
    return $inseq;
	
}


sub canola_b_napus_comparativeGenomic{
	(my $query_file, my $subject_file)=@_;
	#my $query_seq = open_fasta($query_file);
	#concatenate_fasta($query_file,100);
	split_fasta($query_file,100);	
	#my $subject_seq=open_fasta($subject_file);
	split_fasta($subject_file,100)
}

#raffit

#B.Subtillis:  WTGNAANCGNWNNCW  ( N and W stand for any base and A or T respectively) OR 

#WWTGNAARCGNWWWCAWW (where W, R, and N indicate A or T, A or G, and anybase, respectively) 

#C. cellulolyticum: TGTGTACGCGTTTATATT

#===tobin

#In the above sequences, W represents an A or T, N represents any 
#nucleotide (A, C, G or T), and R represents a G or A.  Thus, multiple 
#versions of each sequence may exist within the genome.


#bacillus subtilis
sub degenerate_blast_Ctermitidis{
	my $reg1="WTGNAANCGNWNNCW"; #B.Sub
	my $reg2="WWTGNAARCGNWWWCAWW"; #B.Sub
	my $reg3="TGTGTACGCGTTTATATT";#"RRGAAAANGTTTTCWW"; #C.Cell
	#my @reg=($reg1,$reg2,$reg3);
	#my @reg=($reg1,$reg2);#,$reg3);
	my @reg=($reg3);
	my $deg1="[A|T]"; #W
	my $deg2="[A|C|G|T]"; #N
	my $deg3="[G|A]"; #R
	my @deg=($deg1,$deg2,$deg3);
	my @rep=@reg;
	my $len1=@reg;
	my $len2=@deg;
	for(0..$len1-1){
		my $in1=$_;
		$rep[$in1] =~ s/W/$deg1/g;
		$rep[$in1] =~ s/N/$deg2/g;
		$rep[$in1] =~ s/R/$deg3/g;
	}
	open (DEGFILE, '>C_cell_degenerate_blast.txt');
	print DEGFILE "query\tcontig_id\tstrand\tstart\tend\thit\n";
	my %seq_h=();
	my %rev_h=();
	while (my $seq = $inseq->next_seq){
		#print DEGFILE "=======\ncontig - ".$seq->id."\n=======\n";
		$seq_h{$seq->id}=$seq->seq;
		$rev_h{$seq->id}=$seq->revcom->seq;
	}
	for(@rep){
		my $regex=$_;
		#print DEGFILE "*******\nquery sequence - ".$regex."\n*******\n";
		foreach my $key (keys %seq_h){
			my $cur_seq=$seq_h{$key};
			my $contig_id=$key;
			while($cur_seq =~ /$regex/g){
				my $start=$-[0];
				my $end=$+[0];
				print DEGFILE $regex."\t".$contig_id."\tforward"."\t".$start."\t".$end."\t".(substr $cur_seq, $start,$end-$start+1)."\n"; #-start, +end
				print $regex."\t".$contig_id."\tforward"."\t".$start."\t".$end."\t".(substr $cur_seq, $start,$end-$start+1)."\n"; #-start, +end
			
			}	
			my $rev_seq=$rev_h{$key};
			my $rev_len=length($rev_seq);
			while($rev_seq =~ /$regex/g){
				my $start=$-[0];
				my $end=$+[0];
				my $start_forward=$rev_len-$end+1;
				my $end_forward=$rev_len-$start+1;
				print DEGFILE $regex."\t".$contig_id."\treverse"."\t".$start_forward."\t".$end_forward."\t".(substr $rev_seq, $start,$end-$start+1)."\n"; #-start, +end
				print $regex."\t".$contig_id."\treverse"."\t".$start."\t".$end."\t".(substr $rev_seq, $start,$end-$start+1)."\n"; #-start, +end
			
			}							
		}
	}
	close(DEGFILE);
	
}

#i)  WTGNAANCGNWNNCW
#ii)  WWTGNAARCGNWWWCAWW
#iii)   RRGAAAANGTTTTCWW

#In the above sequences, W represents an A or T, N represents any 
#nucleotide (A, C, G or T), and R represents a G or A.  Thus, multiple 
#versions of each sequence may exist within the genome.

#matching position
#Reference - http://stackoverflow.com/questions/87380/how-can-i-find-the-location-of-a-regex-match-in-perl

sub degenerate_blast{
	
	my $reg1="WTGNAANCGNWNNCW";
	my $reg2="WWTGNAARCGNWWWCAWW";
	my $reg3="RRGAAAANGTTTTCWW";
	my @reg=($reg1,$reg2,$reg3);
	my $deg1="[A|T]"; #W
	my $deg2="[A|C|G|T]"; #N
	my $deg3="[G|A]"; #R
	my @deg=($deg1,$deg2,$deg3);
	my @rep=@reg;
	my $len1=@reg;
	my $len2=@deg;
	for(0..$len1-1){
		my $in1=$_;
		$rep[$in1] =~ s/W/$deg1/g;
		$rep[$in1] =~ s/N/$deg2/g;
		$rep[$in1] =~ s/R/$deg3/g;
	}
	open (DEGFILE, '>degenerate_blast.txt');
	print DEGFILE "query\tcontig_id\tstrand\tstart\tend\thit\n";
	my %seq_h=();
	my %rev_h=();
	while (my $seq = $inseq->next_seq){
		#print DEGFILE "=======\ncontig - ".$seq->id."\n=======\n";
		$seq_h{$seq->id}=$seq->seq;
		$rev_h{$seq->id}=$seq->revcom->seq;
	}
	for(@rep){
		my $regex=$_;
		#print DEGFILE "*******\nquery sequence - ".$regex."\n*******\n";
		foreach my $key (keys %seq_h){
			my $cur_seq=$seq_h{$key};
			my $contig_id=$key;
			while($cur_seq =~ /$regex/g){
				my $start=$-[0];
				my $end=$+[0];
				print DEGFILE $regex."\t".$contig_id."\tforward"."\t".$start."\t".$end."\t".(substr $cur_seq, $start,$end-$start+1)."\n"; #-start, +end
				print $regex."\t".$contig_id."\tforward"."\t".$start."\t".$end."\t".(substr $cur_seq, $start,$end-$start+1)."\n"; #-start, +end
			
			}	
			my $rev_seq=$rev_h{$key};
			my $rev_len=length($rev_seq);
			while($rev_seq =~ /$regex/g){
				my $start=$-[0];
				my $end=$+[0];
				my $start_forward=$rev_len-$end+1;
				my $end_forward=$rev_len-$start+1;
				print DEGFILE $regex."\t".$contig_id."\treverse"."\t".$start_forward."\t".$end_forward."\t".(substr $rev_seq, $start,$end-$start+1)."\n"; #-start, +end
				print $regex."\t".$contig_id."\treverse"."\t".$start."\t".$end."\t".(substr $rev_seq, $start,$end-$start+1)."\n"; #-start, +end
			
			}							
		}
	}
	close(DEGFILE);
}

sub add_header_to_contig{
	open FastaFile,">", $fasta_file."_header.fa";
	#my $header="gnl|UoM|Contig01 [organism=Pseudomonas putida LS46] [strain=LS46] [gcode=11]";
	#my $header="gnl|UoM|Contig1 [organism=Clostridium intestinale strain URNW] [strain=URNW] [country=Canada] [collection-date=2008] [isolation-source=contaminated cellobiose stock bottle] [gcode=11]";
	#my $header="[organism=Caldibacillus debilis] [strain=GB1] [isolation-source=Cellulolytic enrichment culture] [country=Canada] [collection-date=01-Jun-2009]";
	my $header="[organism=Pseudomonas monteilii] [strain=MO2] [isolation source=wastewater from sewage treatment plant] [country=Canada] [collection-date=11-Jul-2012]";
	while (my $seq = $inseq->next_seq) {
		#print FastaFile ">gnl|UoM|".$seq->id." [organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada] [collection-date=Jun-2009] [isolation-source=waste water] [gcode=11]\n";
		print FastaFile ">".$seq->id." ".$header."\n";#[organism=Clostridium intestinale strain URNW] [strain=URNW] [country=Canada] [collection-date=2008] [isolation-source=contaminated cellobiose stock bottle] [gcode=11]\n";
	    print FastaFile $seq->seq."\n";
	}
	close(FastaFile);
}

sub concatenate_fasta{
	(my $input_file,my $gap)=@_;
	open ConcatGapFile,">", $input_file."_Gap".$gap."Ns";
	my $count=1;
	#print ConcatFile ">Clostridium_termitidtis_noGap\n";
	print ConcatGapFile ">".$input_file."_withGap\n";
	my $inseq = open_fasta($input_file);
	while (my $seq = $inseq->next_seq) {
		if ($count>1){
			print ConcatGapFile "N"x$gap;
		}
	    #print ConcatFile $seq->seq;
	    print ConcatGapFile $seq->seq;
	    $count++;
	}
}

sub concatenate_sequence{
	#open ConcatFile,">", $fasta_file."_concatenate_noGap";
	open ConcatGapFile,">", $fasta_file."_concatenate_Gap100Ns";
	my $count=1;
	#print ConcatFile ">Clostridium_termitidtis_noGap\n";
	#print ConcatGapFile ">Clostridium_termitidtis_withGap\n";	
	while (my $seq = $inseq->next_seq) {
		if($count==11){
			#print ConcatGapFile ">99\n".$seq->seq;
			my @codons = $seq->seq =~ /[A-Za-z]{50}/g;
			print ConcatGapFile ">99\n".join("\n",@codons);
		}elsif ($count>11){
			print ConcatGapFile ("N"x50).("N"x50);
			my @codons = $seq->seq =~ /[A-Za-z]{50}/g;
			print ConcatGapFile join("\n",@codons);
		}else{
			print ConcatGapFile ">".$seq->id."\n";
			my @codons = $seq->seq =~ /[A-Za-z]{50}/g;
			print ConcatGapFile join("\n",@codons)."\n";
		}
	    #print ConcatFile $seq->seq;
	    #print ConcatGapFile $seq->seq;
	    $count++;
	}
}
sub grouping_by_id{
	
	open NCBIFile,">", $fasta_file."_length.fa";
	my $base_dir="/local/workspace01/zhangju/p_annotation/Pseudomonas_16s_rRNA";
	my %out_file=();
	while (my $seq = $inseq->next_seq) {
		my $desc=$seq->desc;
		my @arr1=split(/ /,(split(/;/,$desc))[0]);
		my $gName=join('_',@arr1);
		#my $gName="Pseudomonas_".join('_',(split(/ /,(split(/;/,$desc))[0]))[1..2]);
		if(not exists $out_file{$base_dir."/".$gName}){
			open $gName, ">",$base_dir."/".$gName;
			$out_file{$gName}=1;
		}
		print $gName ">".$seq->id." ".$desc."\n";
	    print $gName $seq->seq."\n";
	}
	foreach(keys %out_file){
		close($_);
	}
	
}

sub blast_m8_nohit_gene{
	#$reversed_obj = $seq_obj->revcom;
	my $Scaffold454_file="454Scaffold_concate.fna";
	my $GerhardR_file="CS_Chr_Draft_20121014_2206_R.fasta";
	my $img_scaffold_file="img_all_scaffolds.fna_concatenate_Gap200Ns.fna";
	my $inseq=Bio::SeqIO->new(-file => $Scaffold454_file,
	                               -format => $format
	                              );
	my $Scaffold454_seq=$inseq->next_seq->seq;
	$inseq=Bio::SeqIO->new(-file => $GerhardR_file,
                               -format => $format
                              );
	my $GerhardR_seq=$inseq->next_seq->seq;  
	
    my $seq_obj = Bio::Seq->new(-seq => $GerhardR_seq,-alphabet => 'dna' );
	$reversed_obj = $seq_obj->revcom;
    my $GerhardR_seq_rev=$reversed_obj->seq;
	

	$inseq=Bio::SeqIO->new(-file => $img_scaffold_file,
                               -format => $format
                              );
	my $img_scaffold_seq=$inseq->next_seq->seq;
	
	$seq_obj = Bio::Seq->new(-seq => $img_scaffold_seq,-alphabet => 'dna' );
	$reversed_obj = $seq_obj->revcom;
    my $img_scaffold_seq_rev=$reversed_obj->seq; 
    
	open GeneF, "<img_gene.txt";
	my %gene_h={};
	while(<GeneF>){
		my $line=$_;
		chomp($line);
		my @arr=split("\t",$line);
		$gene_h{$arr[0]}=$arr[1]."\t".$arr[2];
	}
	close(GeneF);
	open hitF, "<gerhard_2470_gene.txt";
	my %hit_h={};
	while(<hitF>){
		my $line=$_;
		chomp($line);
		$hit_h{$line}=1;
	}
	close(hitF);
	my %not_hit_h={};
	open notHitF, ">not_hit_gene.txt";
	foreach(keys %gene_h){
		if(not ( exists $hit_h{$_})){
			print notHitF $_."\n";	
			$not_hit_h{$_}=1;
		}
	}
	close(notHitF);
	my %gene_seq_h={};
	my $img_gene_file = "current_img_allGene.fna_length.fa";
	my $format="Fasta";
	$inseq=Bio::SeqIO->new(-file => $img_gene_file,
	                               -format => $format
	                              );
	open noHitFastaF, ">not_hit_gene.fna";
	open notInNewAssemblyF, ">not_in_newAssembly.fna";
	my $count1=0;
	my $count2=0;
	my $count3=0;
	my $count4=0;
	my $count5=0;
	my $count6=0;
	my $count7=0;
	while (my $seq = $inseq->next_seq) {
		my $seq_id=(split("_",$seq->id))[0];
		my $in_ours=0;
		my $in_new=0;

		if(($GerhardR_seq =~ $seq->seq) || ($GerhardR_seq_rev =~ $seq->seq)){
			#print $seq->id." not in Gerhard R sequence\n";
			$count4++;
			$in_new=1;
		}else{
			$count2++;
			print notInNewAssemblyF ">".$seq_id."\t".$gene_h{$seq_id}."\n";
			print notInNewAssemblyF $seq->seq."\n"
		}

		if(($img_scaffold_seq =~ $seq->seq) || ($img_scaffold_seq_rev =~ $seq->seq)){
			#print $seq->id." not in Gerhard R sequence\n";
			$in_ours=1;
			$count7++;
		}else{
			$count6++;
		}
		if(($in_new==1) && ($in_ours==1)){
			$count5++;
		}
		if(($in_new==0) && ($in_ours==0)){
			print $seq_id."\t";
		}
	
		if (exists $not_hit_h{$seq_id}){
			print noHitFastaF ">".$seq->id."\t".$gene_h{$seq_id}."\n";
			print noHitFastaF $seq->seq."\n";
		}

	}
	close(notInNewAssemblyF);
	close(noHitFastaF);
	#print ("img gene in current assembly - ".$count3."\n");
	#print ("img gene not in current assembly - ".$count1."\n");
	print ("img gene in new assembly - ".$count4."\n");
	print ("img gene not in new assembly - ".$count2."\n");     
	print ("img gene in img scaffold - ".$count7."\n");
	print ("img gene not in img scaffold - ".$count6."\n");                   
	print ("img gene in both assemblies - ".$count5."\n");
	
}
sub blast_m8_img_gene{
	
	open GeneF, "<img_gene.txt";
	my %gene_h={};
	while(<GeneF>){
		my $line=$_;
		chomp($line);
		my @arr=split("\t",$line);
		$gene_h{$arr[0]}=$arr[1]."\t".$arr[2];
	}
	#2514855171
	
	close(GeneF);
	open m8F, "<img_vs_gerhard.m8";
	open m8_geneF, ">img_m8.txt";
	while(<m8F>){
		my $line=$_;
		chomp($line);
		my @arr1=split("\t",$line);
		my @arr2=split("_",$arr1[0]);
		my $product="\t";
		if(exists $gene_h{$arr2[0]}){
			$product=$gene_h{$arr2[0]};
		}
		print m8_geneF $arr2[0]."\t".$arr2[1]."\t".$product."\t".$line."\n";
	}
    close(m8F);
    close(m8_geneF);
		
}
sub sequence_length{
	
	open NCBIFile,">", $fasta_file."_length";
	my $count=0;
	my $countA=0;
	my $countT=0;
	my $countC=0;
	my $countG=0;
	my $countN=0;
	while (my $seq = $inseq->next_seq) {

		print NCBIFile $seq->id."\t".length($seq->seq)."\n";
		$count=$count+length($seq->seq);
		#$countN++ while($seq->seq =~ m/N/g);
		#$countA++ while($seq->seq =~ m/A/g);
		#$countT++ while($seq->seq =~ m/T/g);
		#$countC++ while($seq->seq =~ m/C/g);
		#$countG++ while($seq->seq =~ m/G/g);
		#print $seq->id."\n";
	    #print NCBIFile $seq->seq."\n";

	}
	close(NCBIFile);
	print "total ".$count."\n";
	#print "As ".$countA."\nTs ".$countT."\nCs ".$countC."\nGs ".$countG."\n";
	
}

sub shred_sequence{
	open NCBIFile,">", $fasta_file."_shredded.fa";
	my $strSource="[organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]"; #"[organism=Pseudomonas putida LS46] [strain=LS46][country=CA: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]";
	while (my $seq = $inseq->next_seq) {
		if (length($seq->seq)> 2000){
			my $start=0;
			my $fa_id=$seq->id;
			my $count=1;
			while ($start+300 < length($seq->seq)){
				my $fa_seq = substr($seq->seq, $start, 1500);
				print NCBIFile ">".$seq->id."_".$count."\n";
				print NCBIFile $fa_seq."\n";
				$start=$start+1000;
				$count=$count+1;
			}
		}else{
			print NCBIFile ">".$seq->id."\n";
	    	print NCBIFile $seq->seq."\n";
		}
	}
	close(NCBIFile);
}

sub ncbi_gb1_publication{
	
	open NCBIFile,">", $fasta_file."_NCBI.fsa";
	#my $strSource="[organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]"; #"[organism=Pseudomonas putida LS46] [strain=LS46][country=CA: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]";
	my $strSource="[organism=Caldibacillus debilis] [strain=GB1] [country=Canada] [collection-date=Jul-2009] [isolation-source=marsh]";
	while (my $seq = $inseq->next_seq) {
		my $seq_id=$seq->id;
		if (exists $id_map{$seq->id}){
			$seq_id=$id_map{$seq->id};
		}
		print NCBIFile ">".$seq_id." ".$strSource."\n";
	    print NCBIFile $seq->seq."\n";
	}
	close(NCBIFile);
}

sub ncbi_genome_publication{
	
	open NCBIFile,">", $fasta_file."_NCBI.fsa";
	#my $strSource="[organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]"; #"[organism=Pseudomonas putida LS46] [strain=LS46][country=CA: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]";
	my $strSource="[organism=Caldibacillus debilis] [strain=GB1] [country=Canada] [collection-date=Jul-2009] [isolation-source=marsh]";
	while (my $seq = $inseq->next_seq) {
		my $seq_id=$seq->id;
		if (exists $id_map{$seq->id}){
			$seq_id=$id_map{$seq->id};
		}
		print NCBIFile ">".$seq_id." ".$strSource."\n";
	    print NCBIFile $seq->seq."\n";
	}
	close(NCBIFile);
}

sub ncbi_publication_scaffold{
	

	my %id_map=();
	open IDFile, "</local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/Contig_to_Scaffold.id";
	while(<IDFile>){
		chomp($_);
		my $line=$_;
		my @arr=split(' ',$line);	
		$id_map{$arr[0]}=$arr[1];
	}
	close(IDFile);
	my $format="Fasta";
	my $base_dir = "/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/";
	my $scaffold_old=$base_dir."scaffold_dir/454Scaffolds.fsa_copy";
	my $scaffold_new=$base_dir."scaffold_dir/454Scaffolds.fsa";
    my $scaf_seq=Bio::SeqIO->new(-file => $scaffold_old,
                               -format => $format
                              );
	open ScaffoldFile,">", $scaffold_new;
	my $strSource=" [organism=Pseudomonas putida LS46] [strain=LS46] [country=Canada: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]"; #"[organism=Pseudomonas putida LS46] [strain=LS46][country=CA: Manitoba] [collection-date=Jun-2009] [isolation-source=waste water]";
	while (my $seq = $scaf_seq->next_seq) {
		my $seq_id=$seq->id;
		if (exists $id_map{$seq->id}){
			$seq_id=$id_map{$seq->id};
		}
		print ScaffoldFile ">".$seq_id." ".$strSource."\n";
	    print ScaffoldFile $seq->seq."\n";
	}
	close(NCBIFile);
}

sub ncbi_genome_annotation{
	my $base_dir = "/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/";
	my $ContigFile="/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/fasta_folder/NCBI_publication_allContig.fna_NCBI.fsa";
	open AGPFile, "</local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/454Scaffolds.agp";
	#open FastaFile, "</"
	
	#agp
	
	#scaffold00001   5882767 5907756 301     W       contig00151     1       24990   +
    #scaffold00002   1       5398    1       W       contig00152     1       5398    +
    #scaffold00003   1       3340    1       W       contig00153     1       3340    +
	#####Add each individual contig as an scaffold#####
	
	my %ContigId=();
	my $contigs = Bio::SeqIO->new(-file => $ContigFile,
                               -format => "Fasta"
                              );
	while (my $contig = $contigs->next_seq) {
	    $ContigId{$contig->id}=length($contig->seq);
	}
	
	my %addContig=();
	while(<AGPFile>){
		chomp($_);
		my $line=$_;
		my @arr=split(' ',$line);
		$addContig{$arr[5]}=1;
	}
    close(AGPFile);
	open AGPFile, ">>/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/454Scaffolds.agp";
	open IDMap, ">".$base_dir."Contig_to_Scaffold.id";
	my $start_index=4;
	my %add_scaffold=();
	for my $c_id (keys %ContigId){
		if(not exists $addContig{$c_id}){
			$index = sprintf("%05d", $start_index);
			$scaffold_id="scaffold".$index;
			my $add_line=$scaffold_id."\t1\t".$ContigId{$c_id}."\t1\tW\t".$c_id."\t1\t".$ContigId{$c_id}."\t+\n";
			print AGPFile $add_line;
			print IDMap $c_id."\t".$scaffold_id."\n";
			$add_scaffold{$c_id}=$scaffold_id;
			$start_index=$start_index+1;
		}
	}
	close(IDMap);
	close(AGPFile);
	&ncbi_publication_scaffold();
	open(IMGGeneFile, "</local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/IMG.tbl");
	open NCBIGeneFile,">/local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/scaffold_dir/454Scaffolds.tbl";
	open MiscRNAFile, "</local/workspace01/zhangju/454Data/p_454illuminaPairedf01Unpaired005/assembly/NCBI_publication/miscRNA.description.txt";
	
	my %miscRNA=();
	while(<MiscRNAFile>){
		chomp($_);
		my $line=$_;
		my @arr=split(' ',$line);
		my $size=@arr;
		my $tag=$arr[1];
		$tag =~ s/or/PputLS46_/;
		my $desc=join(' ',@arr[3..$size-1]);
		$miscRNA{$tag}=$desc;
	}
	my %geneH=();
	while(<IMGGeneFile>){
		chomp($_);
		my $line=$_;
		my @arr=split(' ',$line);
		my $size=@arr;
		my $start=$arr[0];
		my $end=$arr[1];
		my $type=$arr[2];
		if ($type eq "miscRNA"){
			$type="ncRNA";
		}
		my $scaffold=$arr[3];
		my $tag=$arr[4];
		$tag =~ s/or/PputLS46_/;
		if ($tag eq "PputLS46_3458"){
			print "check";
		}
		my $product=join(' ',@arr[5..$size-1]);
		if ($product =~ m/(rRNA.+)\)/) {
			my @arr1=split(' ',$1);
			$product = $arr1[1]." ribosomal RNA";
		}
		#tRNA (EC:2.7.7.70) nucleotidyltransferase/poly(A) polymerase (EC:2.7.7.72, EC:3.1.4.-, EC:3.1.3.-)
		my $str_temp=$product;
		my $increase=1;
		
		while($str_temp  =~ m/.+(\(EC:.+\)).+/){
			my $find_str=$1;
			my $result = index($str_temp,$find_str);
			my $pre_str=substr($str_temp,0, $result-1);
			my $pos_str=substr($str_temp, $result+length($find_str)+1,);
			$str_temp=$pre_str." ".$pos_str;
			my @arr=split(', ',substr($1,1,length($1)-2));
			foreach(@arr){
				$geneH{$scaffold}{$tag}{"EC_number"}{$increase}=$_;
				$increase=$increase+1;				
			}
			$product=$str_temp;
		}
		if ($str_temp =~ m/.+(\(EC:.+\))/){
			my $find_str=$1;
			my $result = index($str_temp,$find_str);
			$product=substr($str_temp,0, $result-1);
			my @arr=split(', ',substr($1,1,length($1)-2));
			foreach(@arr){
				$geneH{$scaffold}{$tag}{"EC_number"}{$increase}=$_;
				$increase=$increase+1;				
			}
		}
		#$product=$str_temp;
		if ((substr $start,0,2) eq "F0") {
			$start = (substr $start,2,)
		}
		$geneH{$scaffold}{$tag}{"start"}=$start;
		$geneH{$scaffold}{$tag}{"end"}=$end;
		$geneH{$scaffold}{$tag}{"type"}=$type;
		$geneH{$scaffold}{$tag}{"product"}=$product;
		if ($type eq "miscRNA"){
			$geneH{$scaffold}{$tag}{"product"}=$miscRNA{$tag};
		}
	}
	close(IMGGeneFile);
	for my $scaffold (keys %geneH){
		my $head=$scaffold;
		if(exists $add_scaffold{$scaffold}){
			$head=$add_scaffold{$scaffold};
		}
		print NCBIGeneFile ">Feature ".$head."\n";
		for my $tag (keys $geneH{$scaffold}){
			print NCBIGeneFile $geneH{$scaffold}{$tag}{"start"}."\t".$geneH{$scaffold}{$tag}{"end"}."\tgene\n";
			print NCBIGeneFile "\t\t\tgene\t".$geneH{$scaffold}{$tag}{"product"}."\n";
			print NCBIGeneFile "\t\t\tlocus_tag\t".$tag."\n";
			print NCBIGeneFile $geneH{$scaffold}{$tag}{"start"}."\t".$geneH{$scaffold}{$tag}{"end"}."\t".$geneH{$scaffold}{$tag}{"type"}."\n";
			print NCBIGeneFile "\t\t\tproduct\t".$geneH{$scaffold}{$tag}{"product"}."\n";
			print NCBIGeneFile "\t\t\tprotein_id\tgnl|jgi|".$tag."\n";
			if (exists $geneH{$scaffold}{$tag}{"EC_number"}){
				for my $ec (keys $geneH{$scaffold}{$tag}{"EC_number"}){
					print NCBIGeneFile "\t\t\tEC_number\t".$geneH{$scaffold}{$tag}{"EC_number"}{$ec}."\n";
				}
			}
		}
	}
	close(NCBIGeneFile);
}

sub checking_ncbi_tbl_on_fsa{
	
	my $base_dir="/local/workspace01/zhangju/wc1_annotation/gbk_dir/wc1_contigs";
	our $fsa_file = $base_dir.".fsa";
	our $tbl_file=$base_dir.".tbl";
	my $format="Fasta";
	our $inseq=Bio::SeqIO->new(-file => $fsa_file,
	                               -format => $format
	                              );
	
	my %fsa_seq=();
	while (my $seq = $inseq->next_seq) {
		$fsa_seq{$seq->id}=$seq->seq;
	}
	open FSA, "<", $fsa_file;
	
	while(<FSA>){
		print "avc";
	}
	open tbl_to_fsa,">", $base_dir.".tbl_to_fsa.fa";
	close(tbl_to_fsa);
}

sub break_gap{
	
	open ContigFile, ">".$fasta_file."_contigs.fna";
	my $index=1;
	while (my $seq = $inseq->next_seq) {
		my $str_temp=$seq->seq;
		while($str_temp  =~ m/(N{1,})/){
			my $find_str=$1;
			my $result = index($str_temp,$find_str);
			my $pre_str=substr($str_temp,0, $result);
			my $contigId = "contig".sprintf("%05d", $index);
			#print ContigFile ">".$contigId."\t".length($pre_str)."\n".$pre_str."\n";
			print ContigFile ">".$contigId."\n".$pre_str."\n";
			$str_temp=substr($str_temp,length($pre_str)+length($find_str),);
			$index=$index+1;
		}
		my $contigId = "contig".sprintf("%05d", $index);
		#print ContigFile ">".$contigId."\t".length($str_temp)."\n".$str_temp."\n";
		print ContigFile ">".$contigId."\n".$str_temp."\n";
		$index=$index+1
	}
	close(ContigFile);

}
sub sff_paired_reads{
	
	open R1,">".$fasta_file."_R1.fna";
	open R2,">".$fasta_file."_R2.fna";
	open UnpairedF,">".$fasta_file."_unpaired.fna";
	while (my $seq = $inseq->next_seq) {
		if($seq->id =~ m/left/ ){
			print R1 ">".$seq->id."\n".$seq->seq."\n";
			
		}elsif($seq->id =~ m/right/ ){
			print R2 ">".$seq->id."\n".$seq->seq."\n";
		}else{
			print UnpairedF ">".$seq->id."\n".$seq->seq."\n";
		}
	}
	close(R1);
	close(R2);
	close(UnpairedF);
	
}

sub count_seq{
	my $size=500;
	open OutFile, ">".$fasta_file."_bp".$size.".fna";
	my $index=1;
	while (my $seq = $inseq->next_seq) {

		my $ln=length($seq->seq);
		if($ln>=500){
			print OutFile ">scaffold".$index."_length_".$ln."\n".$seq->seq."\n";
			$index=$index+1;		
		}

	}
	close(OutFile);
}

sub shred{
	open ShredFile, ">".$fasta_file."_shredded.fna";
	my $index=1;
	
	while (my $seq = $inseq->next_seq) {
		my $pos=0;
		my $strtemp=$seq->seq;
		my $ln=length($strtemp);
		if($ln<500){
			next;
		}
		my $start=0;
		while($pos<$ln){
			print ShredFile ">shred".$index."\n";
			print ShredFile substr $strtemp, $pos, 1500;
			print ShredFile "\n";
			$index=$index+1;
			if($pos==0 && $start<26){
				$start=$start+1;
				next;
			}
			$pos=$pos+200;		
		}

	}
	close(ShredFile);
	
}

sub sub_sequence{
	(my $from, my $to)=@_;
	#my $from =shift;
	#my $to=shift;
	#(my $from, my $to)=@ARGV;
	my $seq=$inseq->next_seq->seq;
	my $sub_seq=(substr $seq, $from-1, $to-$from+1);
	print $sub_seq."\n";
}

sub split_fasta{
	
	(my $input_file,my $read_length)=@_;
	my $inseq=open_fasta($input_file);
	open ReadFile, ">".$input_file."_read_bp".$read_length.".fasta";
	while (my $seq = $inseq->next_seq) {
		if (length($seq->seq)> $read_length){
			my $start=0;
			my $fa_id=$seq->id;
			my $count=1;
			while ($start+$read_length < length($seq->seq)){
				my $fa_seq = substr($seq->seq, $start, $read_length);
				print ReadFile ">".$seq->id."_".$count."\n";
				print ReadFile $fa_seq."\n";
				$start=$start+$read_length;
				$count=$count+1;
			}
		}else{
			print ReadFile ">".$seq->id."\n";
	    	print ReadFile $seq->seq."\n";
		}
	}
	close(ReadFile);
	
}

sub split_scaffold{
	
	open contigFile, ">".$fasta_file."_split.contig.fna";
	open AGPFile, ">".$fasta_file."_scaffolds.agp";
	open GapFile,">".$fasta_file."_gap.fna";
	
	my $index=1;
	my $split_count=0;
	my $increase=1;
    my $gap=0;
	while (my $seq = $inseq->next_seq) {
		my $pos=0;
		my $str_temp=$seq->seq;
	    my $total_length=0;
	    my $NCBI_scaffoldId=(split('\|',$seq->id))[0];
	    my $ncbi_scaffold_contig_count=1;
		my $start1=0;
		
	    if(not ($str_temp  =~ m/(N{1,})/)){
	    	$index = sprintf("%05d", $increase);
	    	my $pre_str=$str_temp;
	    	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";
			print AGPFile $agp_sequence;
			print contigFile ">contig".$index."|size".length($pre_str)."\n";
			$index1=0;
			while(($index1+50)<length($pre_str)){
				print contigFile substr($pre_str,$index1,50)."\n";
				$index1=$index1+50;
			}
			print contigFile substr($pre_str,$index1,)."\n";
			$increase++;
			next;
	    }
	    my $adj_base=0;
		while($str_temp  =~ m/(N{1,})/){
			my $find_str=$1;
			#$start1=$total_length;
			$total_length=$total_length+length($1);
			#$total_length=$total_length+$my_gap-1;
			
			my $result = index($str_temp,$find_str);
	
			$start_gap=$result-150;
			if ($start_gap<0){
				$start_gap=0;
			}
			my $len_gap=300;
			if(length($str_temp)-$start_gap+1<300){
				$len_gap=length($str_temp)-$start_gap+1;
			}
			if (length($find_str)<10){
				$gap++;
				print GapFile ">".$NCBI_scaffoldId."_".$gap."\n";
				print GapFile substr($str_temp,$start_gap,$len_gap)."\n";
			}
	
			my $pre_str=substr($str_temp,0, $result);
			$str_temp=substr($str_temp,length($pre_str)+length($find_str),);

			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);
			$index = sprintf("%05d", $increase);
	
			#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";
			my $my_gap=length($find_str);
			if($my_gap<10){
				$my_gap=10;
				$adj_base=$adj_base+10-length($find_str);
				$total_length=$total_length-length($find_str)+10;
				#$start1=$start1+10-length($find_str);
			}else{
				$total_length=$total_length+length($find_str);
			}
			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)+$my_gap)."\t".($ncbi_scaffold_contig_count+1).
			    "\tN\t".$my_gap."\tfragment\tyes\n";
			
			$start1=$start1+length($pre_str)+$my_gap;
			
			print AGPFile $agp_sequence;
			#$total_length=$total_length+1;
			#$start1=$total_length;
			$split_count++;
			print contigFile ">contig".$index."|size".length($pre_str)."\n";
			my $index1=0;
			while(($index1+50)<length($pre_str)){
				print contigFile substr($pre_str,$index1,50)."\n";
				$index1=$index1+50;
			}
			print contigFile substr($pre_str,$index1,)."\n";
	
			$increase=$increase+1;
			$index = sprintf("%05d", $increase);
			$ncbi_scaffold_contig_count=$ncbi_scaffold_contig_count+2;
		}
		print contigFile ">contig".$index."|size".length($str_temp)."\n".$str_temp."\n";
		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;
		$increase=$increase+1;
		$index = sprintf("%05d", $increase);
		
	}
	close(contigFile);
	close(AGPFile);
	close(GapFile);
}
sub img_contig_bowtie_genome{
	open contigFile, ">".$fasta_file."_genome.fna";
	while (my $seq = $inseq->next_seq) {
		print "pause";
		my $str1=(split(':',$seq->desc))[1];
		$str1 =~ s/^\s+//; 
		print contigFile ">".$str1."\n".$seq->seq."\n";
	}
		
	close(contigFile);
	
}

sub get_sequences_byId1{
	
	open ScaffSeq, ">Bnapus_ChromA.fasta";
	open ScaffSeq_C, ">Bnapus_ChromC.fasta";
	while (my $seq = $inseq->next_seq) {
		if ($seq->id =~ "chrA"){
			print ScaffSeq ">".$seq->id."\n".$seq->seq."\n";
		}else{
			print ScaffSeq_C ">".$seq->id."\n".$seq->seq."\n";
		}
	}
	close(ScaffID);
	close(ScaffSeq);
	close(ScaffSeq_C);
}

sub get_sequences_byId_2{
	
	my($id_file)=@_;
	open ScaffID, "<".$id_file;# "<Rsa1.0_00086.1.id";  #"<resistant_gene_locusTag.txt";# "<scaffold_500reads_mismatch012.txt";
	open ScaffSeq,">".$id_file."_fasta.txt";#">Rsa1.0_00086.1.fa";# ">scaffold_500reads_mismatch0to12_sequence.fasta";
	$scaffold_id="";
	while(<ScaffID>){
		chomp($_);
		#my $line=$_;
		my $line = (split / /, $_)[0];
		$scaffold_id=$scaffold_id." ".$line;		
	}
	while (my $seq = $inseq->next_seq) {
		if ($scaffold_id =~ $seq->id){
			print ScaffSeq ">".$seq->id."\n".$seq->seq."\n";
		}
	}

	close(ScaffID);
	close(ScaffSeq);
}

sub get_sequences_byId{
	
	my($id_file,@seq_array)=@_;

	
	open ScaffID, "<".$id_file;# "<Rsa1.0_00086.1.id";  #"<resistant_gene_locusTag.txt";# "<scaffold_500reads_mismatch012.txt";
	open ScaffSeq,">".$id_file."_fasta.txt";#">Rsa1.0_00086.1.fa";# ">scaffold_500reads_mismatch0to12_sequence.fasta";

	$scaffold_id="";
	my %diffGene=();
	while(<ScaffID>){
		chomp($_);
		my $line=$_;
		my @arr1=split("\t",$line);
		$diffGene{$arr1[0]}=join('|',@arr1[0,7,8,9]);
		$scaffold_id=$scaffold_id." ".$line;		
	}
	foreach my $seq ( @seq_array ) {
		if (exists $diffGene{$seq->id}){
			print ScaffSeq ">".$diffGene{$seq->id}."\n".$seq->seq."\n";
		}
	}

=pod	
	$scaffold_id="";
	while(<ScaffID>){
		chomp($_);
		my $line=$_;
		$scaffold_id=$scaffold_id." ".$line;		
	}
	while (my $seq = $inseq->next_seq) {
		if ($scaffold_id =~ $seq->id){
			print ScaffSeq ">".$seq->id."\n".$seq->seq."\n";
		}
	}
=cut
	close(ScaffID);
	close(ScaffSeq);
}

sub tab_merge_fastaLen{
	my %contig=();
	while (my $seq = $inseq->next_seq) {
		$contig{$seq->id}=length($seq->seq);
	}
	
	
	open AGPFile, "<align_mismatch012.reads_3_count";
	open OutFile, ">align_mismatch012.reads_3_count_withScaffoldLength";
	print OutFile "Scaffold\treads\tscaffoldLen\n";
	
	my %addContig=();
	while(<AGPFile>){
		chomp($_);
		my $line=$_;
		my @arr=split(' ',$line);
		print OutFile $arr[0]."\t".$arr[1]."\t".$contig{$arr[0]}."\n";
		$addContig{$arr[0]}=$arr[1];
	}
    close(AGPFile);
    close("AGPFile")
	
}


sub assembly_summary{
	
	my %contig=();
	
	while (my $seq = $inseq->next_seq) {
		$contig{$seq->id}=length($seq->seq);
	}
	my @sort_contig=(sort { $b <=> $a } values(%contig));
	open contigFile, ">scaffold_length.txt";
	while ( my ($key, $value) = each(%contig) ) {
        #print contigFile "$key\t$value\n";
    }
	
	my $total_length=0;
	foreach (@sort_contig)
	{
	    print $_."\t";
	    print contigFile $_."\n";
	    $total_length = $total_length+$_;
	}
	print "\n";
	print "Number of contig - ".@sort_contig."\n";
	print "total length - ".$total_length."\n";
	my $num_contig=@sort_contig;
	print "average length - ".$total_length/$num_contig;
	print "Largest contig - ".$sort_contig[0]."\n";
	my $sum_len=0;
	foreach (@sort_contig)
	{
	    $sum_len = $sum_len+$_;
	    print $_."\n";
	    if($sum_len>$total_length/2){
	    	print "N50 - ".$_."\n";
	    	return;
	    }
	}
	
}

#&checking_ncbi_tbl_on_fsa();
#&blast_m8_img_gene();
#&blast_m8_nohit_gene();
#canola_b_napus_comparativeGenomic("brassicaDB_combinedRapaOleraceas.fa","Brassica_napus_v4.1.chromosomes_genoscope.fa");
#exit 0;


our $fasta_file = shift;
#my $format = shift;
my $format="Fasta";
our $inseq=Bio::SeqIO->new(-file => $fasta_file,
                               -format => $format
                              );
&concatenate_sequence();
#&assembly_summary();
#get_sequences_byId_2("scaffold_id.txt");
exit 0;                       

#&degenerate_blast_Ctermitidis();
#&split_scaffold();
#&sff_paired_reads();
#&ncbi_gb1_publication();
#&tab_merge_fastaLen();
&assembly_summary();
#&shred();
#&add_header_to_contig();
#&img_contig_bowtie_genome();
#&break_gap();

my @seq_array;
while( $seq = $inseq->next_seq() ) {
    push(@seq_array,$seq);
}
&get_sequences_byId("isoform_exp.diff_yes_1vs4",@seq_array);
&get_sequences_byId("isoform_exp.diff_yes_2vs1",@seq_array);
&get_sequences_byId("isoform_exp.diff_yes_3vs1",@seq_array);
&get_sequences_byId("isoform_exp.diff_yes_5vs4",@seq_array);
&get_sequences_byId("isoform_exp.diff_yes_6vs4",@seq_array);

&get_sequences_byId("isoform_exp.diff_yes_7vs1",@seq_array);

&get_sequences_byId("isoform_exp.diff_yes_7vs8",@seq_array);
&get_sequences_byId("isoform_exp.diff_yes_8vs1",@seq_array);
exit 0;                              
my $from =shift;
my $to=shift;

#&sub_sequence($from,$to);#17011,17046);
#print "abc";
#&count_seq();


#&concatenate_sequence();
#&ncbi_genome_annotation();
#&ncbi_genome_annotation();
#&sequence_length();
#&shred_sequence();
#&ncbi_genome_publication()                       
#&concatenate_sequence()                          
# Now that we have a seq stream,
# we need to tell it to give us a $seq.
# We do this using the 'next_seq' method of SeqIO.
