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

use Statistics::Basic qw(:all);
use Bio::SeqIO;

sub bioperl_gff2fasta{
	my $verbose = 0;
 
	## read in gff
	 
	warn "reading GFF\n";
	 
	my %gff;
	 
	open GFF, '<', 'Brapa_gene_v1.5.gff'
	  or die "fail\n";
	 
	while(<GFF>){
	  my ($seqid, undef, undef, $start, $end,
	      undef, undef, undef, $attrs) = split;
	 
	  push @{$gff{$seqid}}, [$start, $end, $attrs];
	}
	 
	warn "OK\n";
	 
	 
	 
	## Do the fasta
	 
	my $seqio = Bio::SeqIO->
	  new( -file => 'Brapa_sequence_v1.5.fa',
	       -format => 'fasta' )
	  or die "double fail\n";
	 
	 
	while(my $sobj = $seqio->next_seq){
	  my $seqid = $sobj->id;
	 
	  unless(defined($gff{$seqid})){
	    warn "no features for $seqid\n";
	    next;
	  }
	 
	  my $seq = $sobj->seq;
	 
	  for(@{$gff{$seqid}}){
	    my ($start, $end, $attrs) = @$_;
	 
	    warn join("\t", $start, $end, $attrs), "\n"
	      if $verbose > 0;
	 
	    my %attrs = split(/=|;/, $attrs);
	 
	    print ">$seqid-". $attrs{"ID"}.
	      "/$start-$end (". ($end-$start+1). ")\n";
	 
	    print substr($seq, $start, $end-$start+1), "\n";
	  }
	 
	  #exit;
	}
	 
	warn "OK\n";
}

sub FPKM_analysis_tab{
	
	my $productFile="product_gff.txt";
	open FILE, $productFile or die $!;
	my %gene_h=();
	my %locus_h=();
	my %translation_h=();
	while(<FILE>){
		my @arr1=split('\t',$_);
		my $gene_id=(split('\.',$arr1[1]))[0];
		my $product_key=0;
		my $locus_key=0;
		my $translation_key=0;
		my $product="";
		foreach my $val(@arr1){
			if($product_key==1){
				$gene_h{$gene_id}=$val;
				$product_key=0;
			}
			if($locus_key==1){
				$locus_h{$gene_id}=$val;
				$locus_key=0;
			}		
			if($translation_key==1){
				$translation_h{$gene_id}=3*(split('\.',$val))[1];
				$translation_key=0;
			}	
			if($val eq "product"){
				$product_key=1;
				next;
			}
			if($val eq "locus_tag"){
				$locus_key=1;
				next;
			}
			if($val eq "translation"){
				$translation_key=1;
				next;
			}
		}
		#$gene_h{$gene_id}=$product;
	}
	close FILE;
	my $pseudoFile="pseudogene_gff.txt";
	open FILE,$pseudoFile or die $!;
	while(<FILE>){
		my @arr1=split('\t',$_);
		my $gene_id=(split('\.',$arr1[1]))[0];
		my $product="pseudogene";
		$gene_h{$gene_id}=$product;
		my $locus_key=0;
		my $translation_key=0;
		foreach my $val(@arr1){
			if($locus_key==1){
				$locus_h{$gene_id}=$val;
				$locus_key=0;
				last;
			}	
			if($translation_key==1){
				$translation_h{$gene_id}=3*(split('\.',$val))[1];
				$translation_key=0;
			}			
			if($val eq "locus_tag"){
				$locus_key=1;
				next;
			}
			if($val eq "translation"){
				$translation_key=1;
				next;
			}
		}
		
	}
	close FILE;
	open OUTFILE, "> locus_tag_product.txt";
	
	foreach $key(sort keys %gene_h){
		print OUTFILE $key."\t".$gene_h{$key}."\n";
	}
	close OUTFILE;
	
	open nonCodingFILE, "nonCoding_gff.txt" or die $!;
	while(<nonCodingFILE>){
		my $line=$_;
		my @arr1=split('\t',$line);
		my $id=(split('=',(split('\.',$arr1[8]))[0]))[1];
		$translation_h{$id}=abs($arr1[3]-$arr1[4])+1;
	}
	close(nonCodingFILE);
	
	open FPKM_SUMFILE, "> FPKM_summary.txt";
	print FPKM_SUMFILE "file\ttotal\tCD\tCD%\trRNA\trRNA%\ttRNA\ttRNA%\n";
	my @files = <FPKM_2013_v2/*fpkm_tracking_v2>;
	my %FPKM_g=();
	foreach $file (@files) {
		open FILE, $file or die $!;
		#FPKM_2013_v2/BM_Rep1_genes.fpkm_tracking_v2
		print FPKM_SUMFILE $file."\t";
		my %fpkm_h=();
		my %line_h=();
		my $head=1;
		my $header="";
		while(<FILE>){
			my $line=$_;
			if($head==1){
				$header=$line;
				$head=0;
				next;
			}
			my $gene_id=(split('\t',$line))[0];
			$line_h{$gene_id}=$line;
			$fpkm_h{$gene_id}=(split('\t',$line))[9];
		}
		close FILE;
		my @names=split('_',(split('\/',$file))[1]);
		my $sample=$names[0].'_'.$names[1];
		#my %fpkm_tmp = (sort {$fpkm_h{$b}<=>$fpkm_h{$a}} keys %fpkm_h);
		my @arr2=values(%fpkm_h);
		$FPKM_g{$sample}=\@arr2;
		
		open OUTFILE, "> ".$file."_withLocusTag.txt";
		#my @headers=split('\t',$header);
		#my $header_size=@headers;
		print OUTFILE "locus_tag\tlength\tproduct\t".$header; #s[0]."\t".join('\t',$headers[1..$header_size-1])."\n";
		my $sum_CD=0;
		my $count_CD=0;
		my $sum_rRNA=0;
		my $count_rRNA=0;
		my $sum_tRNA=0;
		my $count_tRNA=0;
		my $sum_total=0;
		my $count_total=0;
		foreach $key(sort keys %fpkm_h){
			$count_total++;
			$sum_total=$sum_total+$fpkm_h{$key};
			if((substr($gene_h{$key},0,5) eq "tRNA-")&&(length($gene_h{$key})==8)){
				print "tRNA - ".$gene_h{$key}."\t";
				$sum_tRNA=$sum_tRNA+$fpkm_h{$key}/1000*$translation_h{$key};
			}elsif((substr($gene_h{$key},-13) eq "ribosomal RNA")){
				print "rRNA - ".$gene_h{$key}."\t";
				$sum_rRNA=$sum_rRNA+$fpkm_h{$key}/1000*$translation_h{$key};
			}else{
				$sum_CD=$sum_CD+$fpkm_h{$key}/1000*$translation_h{$key};
			}
			print OUTFILE $locus_h{$key}."\t".$translation_h{$key}."\t".$gene_h{$key}."\t".$line_h{$key};
		}
	  	close(OUTFILE);
	  	$sum_total=$sum_CD+$sum_rRNA+$sum_tRNA;
	  	print FPKM_SUMFILE $sum_total."\t".$sum_CD."\t".($sum_CD/$sum_total)."\t".$sum_rRNA."\t".($sum_rRNA/$sum_total)."\t".$sum_tRNA."\t".($sum_tRNA/$sum_total)."\n";
	}
	close(FPKM_SUMFILE);
	open CORFILE, "> FPKM_correlation.txt";
	print CORFILE "Correlation"."\t".join("\t",(keys %FPKM_g))."\n";
	
	foreach $key1(keys %FPKM_g){
		print CORFILE $key1;
		foreach $key2(keys %FPKM_g){
			my $cor1 = correlation($FPKM_g{$key1},$FPKM_g{$key2});
			my @name1=split('_',$key1);
			my @name2=split('_',$key2);
			if($name1[0] eq $name2[0]){
				print CORFILE "\t".(sprintf("%0.4f",$cor1));
			}else{
				print CORFILE "\t"
			}
			
		}
		print CORFILE "\n"
	}
	close(CORFILE);
}

=item FPKM_analysis_gff()

 Arguments:
 	 product_gff.txt - file of "product" lines extracted from annotation gff file
     used for RNA-seq tophat/cufflinks pipeline
 	 pseudo_gff.txt - file of "pseudogene" lines extracted from annotation gff file
     used for RNA-seq tophat/cufflinks pipeline
     FPKM files from RNA-seq analysis
 Outputs: 
 	 locus_tag_product.txt - locus tag with product name
 	 FPKM_summary.txt - FPKM RNA-seq segment ratios of CD, rRNA and tRNA in total
 	 FPKM_withLocusTag.txt - FPKM files adding columns "length of gene"
 	 FPKM_correlation.txt - correlation between replicates of each experient condition
 	 	symbols of correlation table - '=' self 'r' replicates of same condition 'v' exp vs stationary of same condition
 	 

=cut

sub FPKM_analysis_gff{
	
	my $productFile="product_gff.txt";
	open FILE, $productFile or die $!;
	my %gene_h=();
	my %locus_h=();
	my %translation_h=();
	my %non_tag=();
	while(<FILE>){
		my $line1=$_;
		chomp($line1);
		my @arr1=split('\t',$line1);
		my @arr2=split(';',$arr1[8]);
		#for C.Stercorarium
		my $nonTag=False;
		my $nonTagId="";
		my $gene_id="";
		foreach my $val(@arr2){
			if(substr($val,0,10) eq "locus_tag="){
				$gene_id=(split('=',$val))[1];
			}
		}
		#my $gene_id=(split("=",$arr2[2]))[1]; for Putida
		#my $gene_id=(split('\.',(split(';',$arr2[2]))[2]))[0];
		$translation_h{$gene_id}=abs($arr1[3]-$arr1[4])+1;
		my $product="";
		foreach my $val(@arr2){
			if(substr($val,0,8) eq "product="){
				$gene_h{$gene_id}=(split('=',$val))[1];
			}
			if(substr($val,0,10) eq "locus_tag="){
				$locus_h{$gene_id}=(split('=',$val))[1];
			}
		}
		
		if(substr($arr2[0],0,9) ne "ID=PPUTLS"){

			$nonTagId=(split('\.',(split("=",$arr2[0]))[1]))[0];
			if(not exists $non_tag{$nonTagId}){
				$non_tag{$nonTagId}=[$gene_id]
			}else{
				push(@{$non_tag{$nonTagId}},$gene_id);
			}
		}
		
		#$gene_h{$gene_id}=$product;
	}
	close FILE;
#=comment pseudogene
	my $pseudoFile="pseudogene_gff.txt";
	if(-e $pseudoFile){
		open FILE,$pseudoFile or die $!;
		while(<FILE>){
			my $line1=$_;
			chomp($line1);
			my @arr1=split('\t',$line1);
			my $product="pseudogene";
			my @arr2=split(';',$arr1[8]);
			my $gene_id=(split('=',$arr2[1]))[1];
			$gene_h{$gene_id}=$product;
			$translation_h{$gene_id}=abs($arr1[3]-$arr1[4])+1;
			foreach my $val(@arr2){
				if(substr($val,0,8) eq "product="){
					$gene_h{$gene_id}=(split('=',$val))[1];
				}
				if(substr($val,0,10) eq "locus_tag="){
					$locus_h{$gene_id}=(split('=',$val))[1];
				}
			}
		}
		close FILE;
	}
#=cut
	open OUTFILE, "> locus_tag_product.txt";
	
	foreach $key(sort keys %gene_h){
		print OUTFILE $key."\t".$gene_h{$key}."\n";
	}
	close OUTFILE;

	open FPKM_SUMFILE, "> FPKM_summary.txt";
	print FPKM_SUMFILE "file\ttotal\tCD\tCD%\trRNA\trRNA%\ttRNA\ttRNA%\n";
	my @files = <FPKM_2013/*fpkm_tracking>; # for putida, C.termitidis
	#my @files = <FPKM_v2/*fpkm_tracking>; #for C.Stercorarium
	my %FPKM_g=();
	my %line_fpkm=();
	#my $line_sample="locugTag";

	foreach $file (@files) {
		open FILE, $file or die $!;
		#FPKM_2013_v2/BM_Rep1_genes.fpkm_tracking_v2
		print FPKM_SUMFILE $file."\t";
		my %fpkm_h=();
		my %line_h=();

	    my %nonTag_count=();
		my $head=1;
		my $header="";
		while(<FILE>){
			my $line=$_;
			if($head==1){
				$header=$line;
				$head=0;
				next;
			}
			
			my $gene_id=(split('\t',$line))[0];
			if(substr($gene_id,0,6) ne "PPUTLS"){
				if(not exists $nonTag_count{$gene_id}){
					$nonTag_count{$gene_id}=0;
				}else{
					$nonTag_count{$gene_id}=$nonTag_count{$gene_id}+1
				}
				$gene_id=(@{$non_tag{$gene_id}})[$nonTag_count{$gene_id}]
			}
			$line_h{$gene_id}=$line;
			$fpkm_h{$gene_id}=(split('\t',$line))[9];
		}
		close FILE;
		my @names=split('_',(split('\/',$file))[1]);
		#my $sample=$names[0]; putida LS46
		#my $sample=$names[0]."_".$names[1]; #C.stercorarium
		my $sample=$names[0]; #for C.Stercorarium
		#my %fpkm_tmp = (sort {$fpkm_h{$b}<=>$fpkm_h{$a}} keys %fpkm_h);
		my @arr2=values(%fpkm_h);
		
		$FPKM_g{$sample}=\@arr2;
		if($line_sample eq ""){
			$line_sample=(split("_",$sample))[0];
		}else{
			$line_sample=$line_sample."\t".(split("_",$sample))[0];
		}
		
		
		open OUTFILE, "> ".$file."_withLocusTag.txt";
		#my @headers=split('\t',$header);
		#my $header_size=@headers;
		print OUTFILE "locus_tag\tlength\tproduct\t".$header; #s[0]."\t".join('\t',$headers[1..$header_size-1])."\n";
		my $sum_CD=0;
		my $count_CD=0;
		my $sum_rRNA=0;
		my $count_rRNA=0;
		my $sum_tRNA=0;
		my $count_tRNA=0;
		my $sum_total=0;
		my $count_total=0;
		
		foreach $key(sort keys %fpkm_h){
			$count_total++;
			$sum_total=$sum_total+$fpkm_h{$key};
			if((substr($gene_h{$key},0,5) eq "tRNA-")&&(length($gene_h{$key})==8)){
				#print "tRNA - ".$gene_h{$key}."\t";
				$sum_tRNA=$sum_tRNA+$fpkm_h{$key}/1000*$translation_h{$key};
			}elsif((substr($gene_h{$key},-13) eq "ribosomal RNA")){
				#print "rRNA - ".$gene_h{$key}."\t";
				$sum_rRNA=$sum_rRNA+$fpkm_h{$key}/1000*$translation_h{$key};
			}else{
				$sum_CD=$sum_CD+$fpkm_h{$key}/1000*$translation_h{$key};
			}
			#print OUTFILE $locus_h{$key}."\t".$translation_h{$key}."\t".$gene_h{$key}."\t".$line_h{$key};
			my @arr3 = split("\t",$line_h{$key});
			
			if(exists $line_fpkm{$key}){
				$line_fpkm{$key} = $line_fpkm{$key}."\t".$arr3[9];
			}else{
				$line_fpkm{$key} = $arr3[9];
			}
			print OUTFILE $locus_h{$key}."\t".$line_h{$key};
			
			#print $count_total." ";
		}
		print $file."\t".(scalar keys %line_h)."\n";
	  	close(OUTFILE);
	  	$sum_total=$sum_CD+$sum_rRNA+$sum_tRNA;
	  	print FPKM_SUMFILE $sum_total."\t".$sum_CD."\t".($sum_CD/$sum_total)."\t".$sum_rRNA."\t".($sum_rRNA/$sum_total)."\t".$sum_tRNA."\t".($sum_tRNA/$sum_total)."\n";
	}
	close(FPKM_SUMFILE);
	
	open FPKMFILE, "> AllFPKM.txt";
	print FPKMFILE $line_sample."\n";
	foreach $key(sort keys %line_fpkm){
		#print FPKMFILE $key.$line_fpkm{$key}."\n";
		print FPKMFILE $line_fpkm{$key}."\n";
	}
	close(FPKMFILE);
	open CORFILE, "> FPKM_correlation.txt";
	#print CORFILE "Correlation"."\t".join("\t",(keys %FPKM_g))."\n";
	print CORFILE "Correlation";
	foreach $key1(keys %FPKM_g){
		print CORFILE "\t";
		print CORFILE (split('_',$key1))[0];
	}
	print CORFILE "\n";
	foreach $key1(keys %FPKM_g){
		#print CORFILE $key1;
		print CORFILE (split('_',$key1))[0];
		foreach $key2(keys %FPKM_g){
#=comment Putida LS46
			my @name1=split('\.',$key1);
			my @name2=split('\.',$key2);
			my $len1=length($name1[0]);
			my $len2=length($name2[0]);

			my $cor1 = correlation($FPKM_g{$key1},$FPKM_g{$key2});
			if((substr($name1[0],0,$len1-1) eq substr($name2[0],0,$len2-1) && substr($name1[0],-1,1) ne substr($name2[0],-1,1)) && ($name1[1] eq $name2[1])){
				print CORFILE "\tr".(sprintf("%0.4f",$cor1))."r";
			}elsif(($name1[0] eq $name2[0]) && ($name1[1] eq $name2[1])){
				print CORFILE "\t=".(sprintf("%0.4f",$cor1))."=";
			}elsif((substr($name1[0],0,$len1-1) eq substr($name2[0],0,$len2-1)) && ($name1[1] ne $name2[1])){
				print CORFILE "\tv".(sprintf("%0.4f",$cor1))."v";
			}else{
				print CORFILE "\t".(sprintf("%0.4f",$cor1));
				#print CORFILE "\t";
			}
			
		}
		print CORFILE "\n";

=comment Cstercorarum
			my @name1=split('_',$key1);
			my @name2=split('_',$key2);
			my $len1=length($name1[0]);
			my $len2=length($name2[0]);

			my $cor1 = correlation($FPKM_g{$key1},$FPKM_g{$key2});
			
			if(($name1[0] eq $name2[0]) && ($name1[1] ne $name2[1])){
				print CORFILE "\t#".(sprintf("%0.4f",$cor1))."#";
			}elsif(($name1[0] eq $name2[0]) && ($name1[1] eq $name2[1])){
				print CORFILE "\t*".(sprintf("%0.4f",$cor1))."*";
			}else{
				print CORFILE "\t".(sprintf("%0.4f",$cor1));
				#print CORFILE "\t";
			}
			
		}
		print CORFILE "\n";
=cut
	}
	close(CORFILE);
}


sub htseq_count_gff{
	
	my $productFile="product_gff.txt";
	open FILE, $productFile or die $!;
	my %gene_h=();
	my %locus_h=();
	my %translation_h=();
	my %non_tag=();
	while(<FILE>){
		my $line1=$_;
		chomp($line1);
		my @arr1=split('\t',$line1);
		my @arr2=split(';',$arr1[8]);
		#for C.Stercorarium
		my $nonTag=False;
		my $nonTagId="";
		my $gene_id="";
		foreach my $val(@arr2){
			if(substr($val,0,10) eq "locus_tag="){
				$gene_id=(split('=',$val))[1];
			}
		}
		$translation_h{$gene_id}=abs($arr1[3]-$arr1[4])+1;
		my $product="";
		foreach my $val(@arr2){
			if(substr($val,0,8) eq "product="){
				$gene_h{$gene_id}=(split('=',$val))[1];
			}
			if(substr($val,0,10) eq "locus_tag="){
				$locus_h{$gene_id}=(split('=',$val))[1];
			}
		}
		if(substr($arr2[0],0,9) ne "ID=PPUTLS"){

			$nonTagId=(split('\.',(split("=",$arr2[0]))[1]))[0];
			if(not exists $non_tag{$nonTagId}){
				$non_tag{$nonTagId}=[$gene_id]
			}else{
				push(@{$non_tag{$nonTagId}},$gene_id);
			}
		}
		
		#$gene_h{$gene_id}=$product;
	}
	close FILE;
	my $pseudoFile="pseudogene_gff.txt";
	if(-e $pseudoFile){
		open FILE,$pseudoFile or die $!;
		while(<FILE>){
			my $line1=$_;
			chomp($line1);
			my @arr1=split('\t',$line1);
			my $product="pseudogene";
			my @arr2=split(';',$arr1[8]);
			my $gene_id=(split('=',$arr2[1]))[1];
			$gene_h{$gene_id}=$product;
			$translation_h{$gene_id}=abs($arr1[3]-$arr1[4])+1;
			foreach my $val(@arr2){
				if(substr($val,0,8) eq "product="){
					$gene_h{$gene_id}=(split('=',$val))[1];
				}
				if(substr($val,0,10) eq "locus_tag="){
					$locus_h{$gene_id}=(split('=',$val))[1];
				}
			}
		}
		close FILE;
	}
	open OUTFILE, "> locus_tag_product.txt";
	
	foreach $key(sort keys %gene_h){
		print OUTFILE $key."\t".$gene_h{$key}."\n";
	}
	close OUTFILE;

	open FPKM_SUMFILE, "> FPKM_summary.txt";
	print FPKM_SUMFILE "file\ttotal\tCD\tCD%\trRNA\trRNA%\ttRNA\ttRNA%\n";
	my @files = <htseq_count/*count>; # for putida, C.termitidis
	my %FPKM_g=();
	my %line_count=();
    #PPUTLS46_0t25565.r01    8036
    #aat.t01 1371


	foreach $file (@files) {
		open FILE, $file or die $!;
		print FPKM_SUMFILE $file."\t";
		my %count_h=();
		my %line_h=();
	    my %nonTag_count=();

		while(<FILE>){
			my $line=$_;
			if(substr($line,0,10) eq "no_feature"){
				last;
			}

			my $gene_id=(split('\t',((split('\t',$line))[0])))[0];
			if(substr($gene_id,0,6) ne "PPUTLS"){
				if(not exists $nonTag_count{$gene_id}){
					$nonTag_count{$gene_id}=0;
				}else{
					$nonTag_count{$gene_id}=$nonTag_count{$gene_id}+1
				}
				$gene_id=(@{$non_tag{$gene_id}})[$nonTag_count{$gene_id}]
			}
			$line_h{$gene_id}=$line;
			$count_h{$gene_id}=(split('\t',$line))[1];
		}
		close FILE;
		my @names=split('_',(split('\/',$file))[1]);
		my $sample=$names[0]; #for C.Stercorarium
		my @arr2=values(%count_h);
		$FPKM_g{$sample}=\@arr2;
		if($line_sample eq ""){
			$line_sample=(split("_",$sample))[0];
		}else{
			$line_sample=$line_sample."\t".(split("_",$sample))[0];
		}
		open OUTFILE, "> ".$file."_withLocusTag.txt";
		print OUTFILE "locus_tag\tgene_id\tcount\n"; #.$header; #s[0]."\t".join('\t',$headers[1..$header_size-1])."\n";
		foreach $key(sort keys %count_h){

			my @arr3 = split("\t",$line_h{$key});
			
			if(exists $line_count{$key}){
				$line_count{$key} = $line_count{$key}."\t".$count_h{$key};
			}else{
				$line_count{$key} = $arr3[9];
			}
			print OUTFILE $locus_h{$key}."\t".$line_h{$key};
		}
		print $file."\t".(scalar keys %line_h)."\n";
	  	close(OUTFILE);
	}
	close(FPKM_SUMFILE);
	
	open FPKMFILE, "> All_readCount.txt";
	print FPKMFILE $line_sample."\n";
	foreach $key(sort keys %line_count){
		print FPKMFILE $line_count{$key}."\n";
	}
	close(FPKMFILE);
}

&bioperl_gff2fasta();
#&htseq_count_gff();
#&FPKM_analysis_gff();

=for comment
			if(length($key1)!=length($key2)){
				print CORFILE "\t";
				next;
			}
			if($name1[1] ne $name2[1]){
				print CORFILE "\t";
				next;
			}

			if(substr($name1[0],0,$len1-1) ne substr($name2[0],0,$len2-1)){
				print CORFILE "\t";
				next;
			}
			if(substr($name1[0],-1,1) eq substr($name2[0],-1,1)){
				print CORFILE "\t";
				next;				
			}
			#if($name1[0] eq $name2[0]){
			#	print CORFILE "\t";
			#	next;
			#}
=cut
