use Math::Trig;
use Bio::Perl;
use Bio::SeqIO;
use Bio::Annotation::Collection;
use Bio::Annotation::Comment;
use Bio::SeqFeature::Generic;
use Statistics::Descriptive;
use Error qw(:try);
use Number::WithError;
use Math::Symbolic::Custom::ErrorPropagation;
use LWP::Simple;
use Cwd;
use GD::Graph::bars;
use GD::Graph::pie;
use GD::Graph::mixed; 


$dir = getcwd;
system ("rd /s /q pdf");
system ("mkdir pdf");
print "Semi Cleavage Filter for processing results from PeptideProphet\n";
print "developed by the Schilling lab\n";
print "\n THIS SUBVERSION ALSO CREATES A PDF OF MS SPECTRA AS REQUESTED BY MCP \n";
print " \n IT SHOULD BE RUN IN THE ORIGINAL ANALYSIS FOLDER \n SINCE ALL MZXML FILES ARE REQUIRED \n\n";

print "Enter short experiment ID, such as the run numbers.\n";
print "Space not allowed, use_underscore.\n";
print "Example: ST93_102\n";
$experiment = <STDIN>;
chomp $experiment;

print "Ratios will be normalized to treatment/control\n";
print "Control labeled (L)ight or (H)eavy?\n";
$cntrl_label = <STDIN>;
chomp $cntrl_label;
if (($cntrl_label ne "H" )&&($cntrl_label ne "L" ))
{
	print "enter only H or L - script will now terminate\n";
	die;	
}
print "A few words on the type of treament (knock-out, knock-in etc)\n";
$info = <STDIN>;
chomp $info;

# print "Do you wish to control ASAP values with XPRESS?\n (y)es or (n)o?\n ";
# $xpress_cntrl = <STDIN>;
# chomp $xpress_cntrl;

print "Do you wish to create a pdf with MS/MS spectra - caution, takes very long!\n (y)es or (n)o?\n";
$pdf_cntrl = <STDIN>;
chomp $pdf_cntrl;
if ($pdf_cntrl eq "y" )
{
	print "Orbitrap or QSTAR ELITE data?\n O or E?\n";
	$MS_type  = <STDIN>;
	chomp $MS_type ;
}

print "Did you analyze and export with xpress, ppm, hyperscore, precursor neutral mass, and assumed charge (y)es or (n)o?\n";
$asap_ntt = <STDIN>;
chomp $asap_ntt;
if ($asap_ntt ne "y" )
{
	print "Analyze and export with xpress, ppm, and hyperscore - script wil now terminate\n";
	die;	
}

print " Extract semi-specific peptides from a standard  digest (y)es or (n)o?\n";
$filter_ntt = <STDIN>;
chomp $filter_ntt;
if ($filter_ntt eq "y" )
{
	print "Works only if you also exported NTT!\n";	
}
print "Digestion protease (t)rypsin or (g)luc or (c)hymotrypsin? \n";
$digestion_protease = <STDIN>;
chomp $digestion_protease;


print "Did you copy proteome file used for analysis in this directory and named it IPI.fasta (y)es or (n)o?\n";
$asap_ntt = <STDIN>;
chomp $asap_ntt;
if ($asap_ntt ne "y" )
{
	print "copy  proteome file in directory and name it IPI.fasta - script wil now terminate\n";
	die;	
}

open(PROTEINFILE, "IPI.fasta") || die "Cannot open IPI.fasta ";
$PrevFasta="12";
@ipi_array=();
while (<PROTEINFILE>) 
{ 	
	chomp $_;
	$CurrentLine = $_; 
	if($CurrentLine =~m/^>/) 
	{		
		$FASTA = $CurrentLine;
		if (length($PrevFasta) > 2)
		{	
			push (@ipi_array,$PrevFasta);
			push (@ipi_array,$CatSequence);
		}
		$PrevFasta=$FASTA;
		$PrevLine="";
	}
	else
	{	
		chomp($CurrentLine);
		$CurrentLine =~ s/\W//;
		$CurrentLine =~ s/[\s\t\n\r\f\a]//;		
		$CatSequence = $PrevLine.$CurrentLine;
		$PrevLine = $CatSequence;
	}
}
push (@ipi_array,$PrevFasta);
push (@ipi_array,$CatSequence);
close PROTEINFILE;

#START STATISTICS
$tmp=0;
open(interact, "interact.pep.xls") || die "Can not open interact.pep.xls ";
$first_line=<interact>;
$count = 0;
$xpress_pos =0;
$expect_pos = 0;
@first_line_parts=split("\t",$first_line);
foreach(@first_line_parts)
{
	chomp $_;
	if ($_ eq "xpress")
	{
		$xpress_pos = $count;
	}
	if ($_ eq "expect")
	{
		$expect_pos = $count;
	}	
	if ($_ eq "num_tol_term")
	{
		$NTT_pos = $count;
	}	
	if ($_ eq "asap mean")
	{
		$asap_pos = $count;
	}

	$count ++;
}
$stat_xpress = Statistics::Descriptive::Full->new();#for xpress
# $stat_asap = Statistics::Descriptive::Full->new();#for asap
$c1=0;
$c2=0;
while (<interact>) 
{
	chomp $_;
	@line_parts=split("\t",$_);
	$xpress = @line_parts[$xpress_pos];
	$asap=@line_parts[$asap_pos];
	$expect = @line_parts[$expect_pos];
	$xpress = @line_parts[$xpress_pos];	
	$c1++;
	
	# if ($xpress_cntrl eq "n")#remoced fro xpress only
	# {
		$asap=$xpress;
	# }
	
#	if (($xpress > 0) && ($expect < 0.1) && ($xpress ne "[unavailable]") && ($xpress > 0.5) && ($xpress < 2.0) )
	if (($xpress ne "[unavailable]") && ($xpress > 0.33) && ($xpress < 3) && ($asap ne "N/A") && ($asap ne "-1") && ($asap ne "")&& ($asap > 0.33) && ($asap < 3))
	{
	
		$xa=$xpress / $asap;
		if (($xa <= 1.5)&&($xa >= 0.67)) 
		{
			$Fc_xpress=log($xpress)/log(2);  
			$Fc_asap=log($asap)/log(2);  		
			$Fc_xpress =  sprintf ("%.1f", $Fc_xpress);
			$stat_xpress->add_data($Fc_xpress);
			$c2++;
			# $Fc_asap =  sprintf ("%.1f", $Fc_asap);
			# $stat_asap->add_data($Fc_asap);
		}
	}	
}
close interact;

#$mean_Fc_xpress = $stat_xpress->geometric_mean();
$mean_Fc_xpress = $stat_xpress->mode();
$mean_Fc_xpress =  sprintf ("%.2f", $mean_Fc_xpress);
$sd_Fc_xpress = $stat_xpress->standard_deviation();
$sd_Fc_xpress = sprintf ("%.2f", $sd_Fc_xpress);
$xpress_mean = 2**$mean_Fc_xpress;
$xpress_mean = sprintf ("%.2f", $xpress_mean);

#$mean_Fc_asap = $stat_asap->geometric_mean();
# $mean_Fc_asap = $stat_asap->mode();
$mean_Fc_asap =  sprintf ("%.2f", $mean_Fc_xpress);
# $sd_Fc_asap = $stat_asap->standard_deviation();
$sd_Fc_asap = sprintf ("%.2f", $sd_Fc_xpress);
# $asap_mean = 2**$mean_Fc_asap;
$asap_mean = sprintf ("%.2f", $xpress_mean);

$fraction=sprintf ("%.1f", $fraction_tmp);
#END STATISTICS

open(interact, "interact.pep.xls") || die "Can not open interact.pep.xls ";
open(outfile, ">$outname");

$first_line=<interact>;
$count = 0;
$id_pos=0;
$description_pos=0;
$xpress_pos =0;
@first_line_parts=split("\t",$first_line);
foreach(@first_line_parts)
{
	chomp $_;
	if ($_ eq "protein")
	{
		$id_pos = $count;
	}
	if ($_ eq "xpress")
	{
		$xpress_pos = $count;
	}
	if ($_ eq "peptide")
	{
		$peptide_pos = $count;
	}
	if ($_ eq "expect")
	{
		$expect_pos = $count;
	}	
	if ($_ eq "protein_descr")
	{
		$description_pos = $count;
	}
	if ($_ eq "asap mean")
	{
		$asap_pos = $count;
	}
	if ($_ eq "asap error")
	{
		$asap_error_pos = $count;
	}
	if ($_ eq "ppm")
	{
		$ppm_pos = $count;
	}
	if ($_ eq "hyperscore")
	{
		$hyperscore_pos = $count;
	}
	if ($_ eq "probability")
	{
		$pp_pos = $count;
	}
	if ($_ eq "spectrum")
	{
		$spectrum_pos = $count;
	}
	if ($_ eq "assumed_charge")
	{
		$charge_pos = $count;
	}
	if ($_ eq "precursor_neutral_mass")
	{
		$mass_pos = $count;
	}
	if ($_ eq "num_tol_term")
	{
		$ntt_pos = $count;
	}
	
	$count ++;
}

$unique= 0;
$total= 0;

@peptide_info_array1=();
@peptide_info_array2=();
@peptide_final_array=();
@peptide_hypersort_array=();

$condition1=0;	
$condition2=0;	
$condition3=0;	
$condition4=0;	
$condition5=0;	
$condition6=0;	
$condition7=0;

while (<interact>) 
{
	chomp $_;
	$total ++;
	$condition = 0;
	@line_parts=split("\t",$_);
	$id_all = @line_parts[$id_pos];
#	print "$id_all\n";	
	$id_tmp = @line_parts[$id_pos];
	$peptide_raw = @line_parts[$peptide_pos];
#	print "$peptide_raw\n";
	$peptide_clean = $peptide_raw;
	$peptide_clean=~s/\W//g;
	$peptide_clean=~s/\d//g;
	$peptide_clean=~s/[nc]//g;
	$peptide_clean=~s/^[A-Z]//;
	$peptide_clean=~s/[A-Z]$//;
	$pep_length=length($peptide_raw);
	$left = substr($peptide_raw, 0, 1);
	$right = substr($peptide_raw, $pep_length-1, 1);
	$pl2=length($peptide_clean);
	print "\n$total $peptide_clean";
	$description = @line_parts[$description_pos];
	$expect_tmp = @line_parts[$expect_pos];
	$expect = sprintf ("%.4f", $expect_tmp);
	$hyperscore = @line_parts[$hyperscore_pos];
	$ppm = @line_parts[$ppm_pos];	
	$xpress = @line_parts[$xpress_pos];
##################################################################
	$asap = @line_parts[$xpress_pos]; #HERE I SWITCHED TO XPRESS!!!
	$asap_orig=$asap;
	$asap_error = 0.01;#HERE I SWITCHED TO XPRESS!!! WILL NOT BE REPORTED!!!
##################################################################
	$pp = @line_parts[$pp_pos];
	$pp = sprintf ("%.2f", $pp);
	$dta_name = @line_parts[$spectrum_pos];
	$charge = @line_parts[$charge_pos];
	$mass = @line_parts[$mass_pos];
	$mass = sprintf ("%.2f", $mass);
	$ntt = @line_parts[$ntt_pos];

#IF XPRESS_OPTION IS NEGATIVE THEN SET XPRESS EQUALS ASAP
#ALSO DO THIS IN AVERAGING PART
	if ($xpress_cntrl eq "n")
	{
		$asap=$xpress;
	}
	
#determine if peptide is heavy, light, or acetylated
	$N_ter_mod = substr($peptide_raw, 4, 2);
	$pep_type = "";

	$all_modifications = "";
	if ($N_ter_mod eq "29")
	{
		$pep_type = "L";
		$all_modifications = "N-ter +28.03 Da";
	}
	if ($N_ter_mod eq "35")
	{
		$pep_type = "H";
		$all_modifications = "N-ter +34.06 Da";
	}
	if ($N_ter_mod eq "37")
	{
		$pep_type = "H";
		$all_modifications = "N-ter +37.08 Da";
	}
	
	
	
#three types of acetylated peptides: heavy K, light K, no K 
#light K is default	and includes also no K 
	
	if ($N_ter_mod eq "43")
	{
		$pep_type = "Ac";
		$all_modifications = "N-ter + 42.01 Da";
	}
	
	if (($N_ter_mod eq "43") && (index($peptide_raw, "162.16") > -1))
	{
		$pep_type = "Ac_heavy";
		$all_modifications = "N-ter + 42.01 Da";
	}

	if (index($peptide_raw, "162.16") > -1)
	{
		$all_modifications = $all_modifications . ", K +34.06 Da";
	}

	if (index($peptide_raw, "156.13") > -1)
	{
		$all_modifications = $all_modifications . ", K +28.03 Da";
	}

	if (index($peptide_raw, "164.17") > -1)
	{
		$all_modifications = $all_modifications . ", K +37.08 Da";
	}
	

	if (index($peptide_raw, "160.03") > -1)
	{
		$all_modifications = $all_modifications . ", C +57.02 Da";
	}

	
#strategy: a number of conditions must be fulfilled: just sum it up and only print if that number has been reached
	print " ";
#condition 1: expect-value < 2.0
	if ($expect < 2.0)
	{
		$condition ++;
		print "1";
		$condition1 ++;
	}
	
#condition 2: C-terminal residue is tryptic
	$c_ter_tryptic = 0;
	$last = substr $peptide_clean,-1,1;
	if ((($last eq "K")&&($filter_ntt eq "n" ) && ($digestion_protease eq "t")) || (($last eq "R")&&($filter_ntt eq "n" )&& ($digestion_protease eq "t")))
	{
		$condition ++;
		print "2";
		$condition2 ++;
		$c_ter_tryptic = 1;
	}
	
	if ((($last eq "E")&&($filter_ntt eq "n" ) && ($digestion_protease eq "g")))
	{
		$condition ++;
		print "2";
		$condition2 ++;
		$c_ter_tryptic = 1;
	}


#exception for semi tryptic analysis
	$subsequent_residue=substr($peptide_raw,length($peptide_raw)-1,1); 	
	
	$semi_tryp_exception=0;
	print " $subsequent_residue ";
	

	if (($filter_ntt eq "y" ) && ($ntt eq 1) && ($subsequent_residue ne "-")) #also reject semi-spec peptides from protein N and C-termini
	{
		$condition ++;
		print "2e";
		$condition2 ++;
		$semi_tryp_exception=1;
	}

#BUILD A TRAP here recognize absence of k and non determined a and x	
	$exception =0;
	if (($xpress eq "[unavailable]") && ($peptide_clean !=~ m/"K"/))
	{
		print "E";
		$exception=1;	
	}	
	
#cond 3 asap exists	
	if (($asap ne "N/A") && ($asap ne "-1") && ($asap != "-1") ) #EXCEPTION NOT NEEDED
	{
		$condition ++;
		print "3";
		$condition3 ++;
	}
#make this optional: $xpress_option eq "y" not yet implemented
#condition 4: xpress exists 
	if (($xpress ne "N/A") && ($xpress ne "-1") && ($xpress != "-1")) #XPRESS EXISTS
	{
		$condition ++;
		print "4";
		$condition4 ++;
	}
	$xpress = $xpress / $xpress_mean;
	$asap = $asap / $asap_mean; 


#make this optional: $xpress_option eq "y" maybe just accept 
#	either both are over 2.5 or below 0.4, x/a should not be larger than 2 or smaller than 0.5

	if ($exception == 0)
	{
		if ($xpress <= 0.001)
		{
			$xpress = 0.001;
		}	
		if ($asap <= 0.001)
		{
			$asap = 0.001;
		}
		if ($asap_orig <= 0.001)
		{
			$asap_orig = 0.001;
		}
		
		$error_ratio = $asap_error / $asap;
		$xa=$xpress / $asap;
		if ((($xpress >= 2.5)&&($asap >= 2.5)) || (($xpress <= 0.4)&&($asap <= 0.4)) || (($xa <= 2)&&($xa >= 0.5))) 
		{
			$condition ++;
			print "5";
			$condition5 ++;
		}
		if ($cntrl_label ne "H" )
		{
			$xpress = 1  / $xpress;
			$asap = 1 / $asap;
			$asap_error = $asap * $error_ratio;	
			$asap_orig = 1 / $asap_orig;
		}	
		if ($xpress <= 0.001)
		{
			$xpress = 0.001;
		}	
		if ($asap <= 0.001)
		{
			$asap = 0.001;
		}
		if ($asap_orig <= 0.001)
		{
			$asap_orig = 0.001;
		}		
		$xpress = sprintf ("%.3f", $xpress);	
		$asap = sprintf ("%.3f", $asap);
		$asap_orig = sprintf ("%.3f", $asap_orig);
		$asap_error = sprintf ("%.3f", $asap_error);
	}
	if ($exception == 1)
	{
		$condition ++;
	}
	
#condition 6
#p1prime is NOT a proline  - should also apply to acetylation
	$p1prime = substr ($peptide_clean,0,1);
	if ($p1prime ne "P")
	{
		$condition ++;
		print "6";
		$condition6 ++;
	}

	# #condition 7 - only when used for semi-tryptic filter
	# if ($filter_ntt eq "y" )	
	# {
		# $condition = $condition - 1;
		# print "7";
		# $condition7 ++;
		# if ($ntt eq "1" )
		# {
			# $condition ++;
		# }
	# }

	print " L520 condition-number is $condition";
			
	if ($condition eq 6)	
	{	
		if ($exception == 1)
		{
			$asap = "n.d.";
			$asap_orig = "n.d.";
			$xpress = "n.d.";
		}
		print "ntt_status = $semi_tryp_exception\n";
		$info_entry = $peptide_clean . "__" . $xpress . "__" . $asap . "__" . $id_all . "__" . $pp . "__" . $hyperscore . "__" . $ppm . "__" . $asap_error . "__" . $peptide_raw . "__" . $dta_name . "__" . $pep_type . "__" .$charge . "__" . $mass . "__" . $all_modifications . "__" . $asap_orig . "__" . $semi_tryp_exception;
		push (@peptide_info_array1,$info_entry);
		push (@peptide_info_array2,$info_entry);	
	}
}


$counter2 =0;
foreach $entry1 (@peptide_info_array1) 
{
	$counter2 ++;
}
print "\n $counter2 peptides passed first level\n";
print "\n c1: $condition1, c2: $condition2, c3: $condition3, c4: $condition4, c5: $condition5 c6: $condition6 c7: $condition7\n";


foreach $entry1 (@peptide_info_array1) 
{
	@entry_parts1=split("__",$entry1);
	$array_peptide = @entry_parts1[0];
	$array_xpress=@entry_parts1[1];
	$array_asap=@entry_parts1[2];
	$array_ids=@entry_parts1[3];
	$array_pp=@entry_parts1[4];
	$array_hyperscore = @entry_parts1[5];	
	$array_ppm=@entry_parts1[6];
	$array_asap_error=@entry_parts1[7];
	$repeated_peptide = 0;
	$repeated_peptide_highest_score = 1;
	
	foreach $entry2 (@peptide_info_array2) 
	{
		@entry_parts2=split("__",$entry2);
		$array_peptide2 = @entry_parts2[0];
		$array_hyperscore2 = @entry_parts2[5];
		if($array_peptide eq $array_peptide2)
		{	
			$repeated_peptide ++;	

			if ($array_hyperscore2 > $array_hyperscore)
			{	
				$repeated_peptide_highest_score = 0;
			}
		}	
	}
	if($repeated_peptide ==1) #must be found at least once in second array...
	{
		push (@peptide_final_array,$entry1);
	}
	if(($repeated_peptide > 1)&&($repeated_peptide_highest_score == 1))
	{
#peptides with identical hyperscore would show up twice - this must be avoided
#so check that peptide is not already present
		$peptide_already_in_list = 0;
		foreach $entry3 (@peptide_final_array) 
		{
			@entry_parts3=split("__",$entry3);
			$array_peptide3 = @entry_parts3[0];
			if($array_peptide eq $array_peptide3)
			{
				$peptide_already_in_list = 1;
			}
		}
		if($peptide_already_in_list == 0)
		{
			push (@peptide_final_array,$entry1);
		}
	}
}

$stat_percentiles = Statistics::Descriptive::Full->new();#for percentiles
open(msms, ">msms.txt");
foreach $entry (@peptide_final_array) 
#$asap_error7 . "__" . $peptide_raw8 . "__" . $dta_name9 . "__" . $pep_type10;
{
	@entry_parts=split("__",$entry); 
	$array_peptide = @entry_parts[0];
	$array_xpress = @entry_parts[1];
	$array_asap = @entry_parts[2];
	$id_all = @entry_parts[3];
	$array_pp =@entry_parts[4];
	$array_hyperscore = @entry_parts[5];
	$array_ppm = @entry_parts[6];
	$array_asap_error = @entry_parts[7];
	$array_dta = @entry_parts[9];
	$array_pep_type = @entry_parts[10];
	$array_charge = @entry_parts[11];
	$array_mass = @entry_parts[12];	
	$array_modifications = @entry_parts[13];
	$array_asap_orig = @entry_parts[14];
	$array_ntt_status = @entry_parts[15];
	
	@dta_parts=split("-s1",$array_dta); 
	$dta_first_half = @dta_parts[0] . "-s1";
	
	if ($MS_type eq "O")
	{
		@dta_parts=split(/\./,$array_dta); 
		$dta_first_half = @dta_parts[0];
	}
	
#make spectrum pdfs
#still need to determine every modified aa and its position!
	
	$residue = 'C';
	$C_mods ="";
	$offset = 0;
	$result = index($array_peptide, $residue, $offset);
	while ($result != -1)
	{
		print "Found $residue at $result\n";
		$offset = $result + 1;
		$result = index($array_peptide, $residue, $offset);
		$C_mods = $C_mods ."&Mod" . $offset . "=160.0306";
	}
	
	if($array_pep_type eq "H")
	{
		$residue = 'K';
		$K_mods ="";
		$offset = 0;
		$result = index($array_peptide, $residue, $offset);
		while ($result != -1)
		{
			print "Found $residue at $result\n";
			$offset = $result + 1;
			$result = index($array_peptide, $residue, $offset);
			$K_mods = $K_mods ."&Mod" . $offset . "=162.1581";
		}
		$outfile = "./pdf/" . $array_peptide . "_heavy". ".pdf";
#		$input = "http://localhost/tpp-bin/plot-msms.cgi?MassType=1&NumAxis=1&ShowB=1&ShowY=1&ShowY2=1&ModN=35.0709&Mod12=162.1581&Pep=" . $array_peptide .  "&Dta=" . $dir . "/" . $dta_first_half .  "/" . $array_dta .  ".dta&Global_Mod=CS57.021500000000003&Global_Mod=KS34.063099999999999&Global_Mod=CV-17.026499999999999&Global_Mod=EV-18.0106&Global_Mod=QV-17.026499999999999";
		$input = "http://localhost/tpp-bin/plot-msms.cgi?MassType=1&NumAxis=1&ShowB=1&ShowY=1&ShowY2=1&ModN=35.0709" .$K_mods. $C_mods."&Pep=" . $array_peptide .  "&Dta=" . $dir . "/" . $dta_first_half .  "/" . $array_dta .  ".dta&Global_Mod=CS57.021500000000003&Global_Mod=KS34.063099999999999&Global_Mod=CV-17.026499999999999&Global_Mod=EV-18.0106&Global_Mod=QV-17.026499999999999";
	}
	
	if($array_pep_type eq "L")
	{
		$residue = 'K';
		$K_mods ="";
		$offset = 0;
		$result = index($array_peptide, $residue, $offset);
		while ($result != -1)
		{
			print "Found $residue at $result\n";
			$offset = $result + 1;
			$result = index($array_peptide, $residue, $offset);
			$K_mods = $K_mods ."&Mod" . $offset . "=156.1263";
		}
		$outfile = "./pdf/" . $array_peptide . "_light". ".pdf";
		$input = "http://localhost/tpp-bin/plot-msms.cgi?MassType=1&NumAxis=1&ShowB=1&ShowY=1&ShowY2=1&ModN=29.0391" .$K_mods. $C_mods. "&Pep=" . $array_peptide .  "&Dta=" . $dir . "/" . $dta_first_half .  "/" . $array_dta .  ".dta&Global_Mod=CS57.021500000000003&Global_Mod=KS34.063099999999999&Global_Mod=CV-17.026499999999999&Global_Mod=EV-18.0106&Global_Mod=QV-17.026499999999999";
	}


	if($array_pep_type eq "Ac")
	{
		$residue = 'K';
		$K_mods ="";
		$offset = 0;
		$result = index($array_peptide, $residue, $offset);
		while ($result != -1)
		{
			print "Found $residue at $result\n";
			$offset = $result + 1;
			$result = index($array_peptide, $residue, $offset);
			$K_mods = $K_mods ."&Mod" . $offset . "=156.1263";
		}
		$outfile = "./pdf/" . $array_peptide . "_acetylated". ".pdf";
		$input = "http://localhost/tpp-bin/plot-msms.cgi?MassType=1&NumAxis=1&ShowB=1&ShowY=1&ShowY2=1&ModN=43.0184" .$K_mods. $C_mods. "&Pep=" . $array_peptide .  "&Dta=" . $dir . "/" . $dta_first_half .  "/" . $array_dta .  ".dta&Global_Mod=CS57.021500000000003&Global_Mod=KS34.063099999999999&Global_Mod=CV-17.026499999999999&Global_Mod=EV-18.0106&Global_Mod=QV-17.026499999999999";
	}
	
	if($array_pep_type eq "Ac_heavy")
	{
		$residue = 'K';
		$K_mods ="";
		$offset = 0;
		$result = index($array_peptide, $residue, $offset);
		while ($result != -1)
		{
			print "Found $residue at $result\n";
			$offset = $result + 1;
			$result = index($array_peptide, $residue, $offset);
			$K_mods = $K_mods ."&Mod" . $offset . "=162.1581";
		}
		$outfile = "./pdf/" . $array_peptide . "_heavy_acetylated". ".pdf";
		$input = "http://localhost/tpp-bin/plot-msms.cgi?MassType=1&NumAxis=1&ShowB=1&ShowY=1&ShowY2=1&ModN=43.0184" .$K_mods. $C_mods. "&Pep=" . $array_peptide .  "&Dta=" . $dir . "/" . $dta_first_half .  "/" . $array_dta .  ".dta&Global_Mod=CS57.021500000000003&Global_Mod=KS34.063099999999999&Global_Mod=CV-17.026499999999999&Global_Mod=EV-18.0106&Global_Mod=QV-17.026499999999999";
	}

	print msms "$array_pep_type array_dta is $array_dta dta_first_half is $dta_first_half $input \n";
	$command = "wkhtmltopdf " . $input . " " . $outfile;

	if ($pdf_cntrl eq "y")
	{	
		system($command);
	}
	
	if (($array_xpress ne "n.d.")&&($array_asap ne "n.d."))
	{		
		$Fc_asap=log($array_asap)/log(2);
		$Fc_asap_orig=log($array_asap_orig)/log(2);
		$Fc_asap = sprintf ("%.3f", $Fc_asap);
		$stat_percentiles->add_data($Fc_asap);
		$error_ratio = $array_asap_error / $array_asap;
		if ($Fc_asap gt 0)
		{
			$Fc_asap_error=$Fc_asap * $error_ratio;
		}	
		if ($Fc_asap lt 0)
		{
			$Fc_asap_error=$Fc_asap * $error_ratio;
			$Fc_asap_error=$Fc_asap_error*-1;
		}
		if ($Fc_asap_error lt 0.001)
		{
			$Fc_asap_error =  0.001;
		}	
		$Fc_asap_error = sprintf ("%.3f", $Fc_asap_error);
	}
	if (($array_xpress eq "n.d.")&&($array_asap eq "n.d."))
	{
		$Fc_xpress="n.d.";
		$Fc_asap = "n.d.";
		$error_ratio = "n.d.";
		$Fc_asap_error="n.d.";
		$Fc_asap_orig="n.d.";
	}
#begin PICS routine
	$PrimeRes = $array_peptide; #keep use PrimeRes variable in PICS routine
	chomp($array_peptide);
	if($array_peptide =~m/^[ACDEFGHIKLMNPQRSTVWY]/) 
	{
		foreach(@ipi_array)
		{
			chomp $_;
			$CurrentLine=$_;
			$seq_matches = 0;
			$PrimeLength=length($array_peptide);
			if($CurrentLine =~ m/^[ACDEFGHIKLMNPQRSTVWY]/) 
			{		
				$CurrentLine =~ s/\W//;
				while($CurrentLine =~ m/($PrimeRes)/g)
				{
					$POSITION=pos($CurrentLine);
					$POSITION_REPORT=$POSITION-$PrimeLength -1;
					$ProteinLength = length($CurrentLine);
					if (($POSITION - $PrimeLength) < 10)
					{
						$START = 0;
						$SubstrLength = $POSITION - $PrimeLength;
					}
					else
					{
						$START= $POSITION - $PrimeLength -10;
						$SubstrLength = 10;						
					}
					$STARTfollowing=$POSITION;
					$Preceeding10=substr($CurrentLine,$START,$SubstrLength);
					if ($array_ntt_status == 1 )
					{
						$Following10=substr($CurrentLine,$STARTfollowing,10);
					}
					else #for real TAILS experiments
					{
						$Following10="";
					}
					
				}
			}
		}
	}	
#end PICS routine
#		print outfile " \t \n$experiment \t$NonPrime\t$array_peptide\t$array_asap\t$array_asap_error\t$Fc_asap\t$Fc_asap_error\t$array_xpress\t$Fc_xpress\t$array_hyperscore\t$array_ppm\t$P1prime\t$ProteinLength";

#jan17
		# if ($array_ntt_status == 1 )
		# {
			# $NonPrime = "tryptic N-terminus ";
		# }
		
		
		$digestion_protease; #can be t or g
		#i need last character of array peptide
		# and an exception for true c-termini - but that is second priority
		
		$array_peptide_length=length($array_peptide);
		$preceeding10_length=length($Preceeding10);
		$last_residue_peptide=substr($array_peptide,$array_peptide_length-1,1);
		$last_residue_preceeding10=substr($Preceeding10,$preceeding10_length-1,1);

		$MS="";
		$DB="";
		$accept="no";
		if (($digestion_protease eq "c") && ($last_residue_peptide eq "F")&& ($last_residue_preceeding10 ne "M"))
		{
			$MS=prime;
			$DB=nonprime;
			$accept = yes;
		}
		if (($digestion_protease eq "c") && ($last_residue_peptide eq "W")&& ($last_residue_preceeding10 ne "M"))
		{
			$MS=prime;
			$DB=nonprime;
			$accept = yes;
		}
		if (($digestion_protease eq "c") && ($last_residue_peptide eq "Y")&& ($last_residue_preceeding10 ne "M"))
		{
			$MS=prime;
			$DB=nonprime;
			$accept = yes;
		}
		if (($digestion_protease eq "c") && ($last_residue_peptide eq "L")&& ($last_residue_preceeding10 ne "M"))
		{
			$MS=prime;
			$DB=nonprime;
			$accept = yes;
		}
		
		if (($digestion_protease eq "c") && ($last_residue_preceeding10 eq "F")&& ($last_residue_peptidene ne "M"))
		{
			$MS=nonprime;
			$DB=prime;
			$accept = yes;
		}		
		if (($digestion_protease eq "c") && ($last_residue_preceeding10 eq "W")&& ($last_residue_peptidene ne "M"))
		{
			$MS=nonprime;
			$DB=prime;
			$accept = yes;
		}
		if (($digestion_protease eq "c") && ($last_residue_preceeding10 eq "L")&& ($last_residue_peptidene ne "M"))
		{
			$MS=nonprime;
			$DB=prime;
			$accept = yes;
		}
		if (($digestion_protease eq "c") && ($last_residue_preceeding10 eq "Y")&& ($last_residue_peptidene ne "M"))
		{
			$MS=nonprime;
			$DB=prime;
			$accept = yes;
		}
		
		if (($digestion_protease eq "g") && ($last_residue_peptide eq "E") && ($last_residue_preceeding10 ne "D") )
		{
			$MS=prime;
			$DB=nonprime;
			$accept = yes;
		}
		if (($digestion_protease eq "g") && ($last_residue_preceeding10 eq "E") && ($last_residue_peptide ne "D"))
		{
			$MS=nonprime;
			$DB=prime;
			$accept = yes;
		}
		
		if (($digestion_protease eq "t") && (($last_residue_peptide eq "K")||($last_residue_peptide eq "R")) && ($last_residue_preceeding10 ne "K")&& ($last_residue_preceeding10 ne "R"))
		{
			$MS=prime;
			$DB=nonprime;
			$accept = yes;
		}
		if (($digestion_protease eq "t") && (($last_residue_preceeding10 eq "K")||($last_residue_preceeding10 eq "R")&&($last_residue_peptide ne "R"))&&($last_residue_peptide ne "K"))
		{
			$MS=nonprime;
			$DB=prime;
			$accept = yes;		
		}

		if ($filter_ntt eq "n" )		
		{
			$accept = yes;		
		}
		
#		$ALL_POSITIONS=99;
		print "accept = $accept $Preceeding10-$array_peptide-$Following10 $last_residue_preceeding10 $last_residue_peptide $MS $DB(XPRESS=$array_xpress) passed filters. \nPossible proteins are:\n";
		if ($accept eq "yes")
		{
			print("\n accepted for output \n");
			$entry_for_averaging_array= $array_peptide. "x_x" . $Preceeding10 . "x_x" . $Fc_asap . "x_x" . $Fc_asap_error . "x_x" . $array_hyperscore . "x_x" .$array_ppm . "x_x" . $array_charge . "x_x" . $array_mass . "x_x" .	$array_modifications  . "x_x" .$Fc_asap_orig . "x_x" .$POSITION_REPORT . "x_x" .$ProteinLength. "x_x" . $Following10 . "x_x" . $MS . "x_x" . $DB;
			@swiss_tmp=();
			$tmp_print="";
			$tmp_info="";
			@individual_ids=split(",",$id_all);
			foreach(@individual_ids)
			{	
				$_=~s/\-/x/g;
				$_=~s/\W/,/g;
				$_=~s/x/-/g;			
				@id_split = split(",",$_); 
				$tmp_print=$tmp_print."\t".@id_split[1];
				$tmp_info= $tmp_info . @id_split[1]."x_x";
				print "@id_split[1] \n";	
			} 
	#		print outfile "$tmp_print";
			$entry_for_averaging_array= $entry_for_averaging_array . "x_x" . $tmp_info;
			push (@array_for_averaging, $entry_for_averaging_array);
			$entry_for_averaging_array="";
			$tmp_info="";
		}	
#	}
	print "\n";	
} 
#print outfile " \t \n";
close interact;
#close outfile;


$low_percentile_cutoff=$stat_percentiles->percentile(20);
$high_percentile_cutoff=$stat_percentiles->percentile(80);

@array_for_averaging=sort(@array_for_averaging);
$prev_prime="";
$outname_non_red = $experiment."_"."TAILS_peptides_non-redundant.xls";
open(outfile_non_red, ">$outname_non_red");
#print outfile_non_red "Mode Fc-value for XPRESS = $mean_Fc_xpress \n This equals an XPRESS value of $xpress_mean\n";
print outfile_non_red "Mode Fc-value for XPRESS = $mean_Fc_asap \n This equals an ASAP value of $asap_mean\n";
print outfile_non_red "XPRESS  normalized to (knock-out, knock-down, treated, overexpressed, etc) over control.\n";
print outfile_non_red "Control was stated as being: $cntrl_label\n";
print outfile_non_red "Infos on the type of treatment and control: $info.\n";
# print outfile_non_red "Fc value is based on XPRESS. ASAPratio was verified by XPRESS (yes/no):$xpress_cntrl\n";
print outfile_non_red " Quantiles: 20 percent lowest Fc values are 'less'; 20 percent highest Fc values are 'more'\n";
print outfile_non_red "For information asap without normalization over normal distribution is also shown. This value is however normalized to treated over control.\n";
print outfile_non_red "Experiment\tMS\tDB\tnonprime\tprime\tId1\tId2\tId3\tname\tFc normalized\tFc error\tquantile\tFc pre-normalized\tPosition\tProtein length\tModifications\tHyperscore\tppm\tcharge\tmass\tAlternate names\tGO compartment\tGO process\tGO function\tSP localization\tsignal peptide\tsecreted or cell surface \n";

#for the histogram
$h1=0;
$h2=0;
$h3=0;
$h4=0;
$h5=0;
$h6=0;
$h7=0;
$h8=0;
$h9=0;
$total_ids=0;
$sec_ids=0;
foreach(@array_for_averaging)
{
	@linesplit=split("x_x",$_);
	$prime= @linesplit[0]; 
	$nonprime= @linesplit[1]; 
	$Fc = @linesplit[2];
	$Fc_error = @linesplit[3];
	$hyperscore= @linesplit[4];
	$ppm= @linesplit[5];
	$charge = @linesplit[6];
	$mass = @linesplit[7];
	$modifications = @linesplit[8];
	$Fc_orig = @linesplit[9];
	$Position = @linesplit[10];
	$Prot_Length = @linesplit[11];	
	$Following10 = @linesplit[12];
	$MS= @linesplit[13];
	$DB= @linesplit[14];	
	$id1=@linesplit[15];
	$id2=@linesplit[16];
	$id3=@linesplit[17];
	
	if($MS eq "nonprime")
	{	
		$nonprime=$prime;
		$prime=$Following10;
	}

	$id = $id1 . "\t" . $id2 . "\t" . $id3;
	$AltNames="";
	$description="";
	$GO_all_C = "";
	$GO_all_P = "";
	$GO_all_F = "";
	$output_loc = "";
	$SubName="";
	$signal = " ";
	$gene_name = " ";
	$prints_family = " ";
	$pfam= " ";
	$interpro= " ";

	$url="http://www.uniprot.org/uniprot/" .$id1. ".txt";
	$page = get($url);
	print "$url \n";
	open(tmp, ">tmp.txt");
	print tmp $page;
	close tmp;
	open(tmp2, "tmp.txt");
	while (<tmp2>)
	{
		chomp $_;
		if($_ =~m/^DE   SubName: Full=/)
		{	
			@subname_split=split("Full=",$_);
			$SubName = @recname_split[1];
		}
		if($_ =~m/^DE   RecName: Full=/)
		{	
			@recname_split=split("Full=",$_);
			$description = @recname_split[1];
		}			
		if($_ =~m/^DE   AltName: Full=/)
		{	
			@altname_split=split("=",$_);
			$AltNames = $AltNames . @altname_split[1];
		}
		if($_ =~m/^DR   GO;/)
		{	
			@go_split_level1=split(";",$_);
			@go_split_level2=split(":",@go_split_level1[2]);# 0= category (C, P, F), 1 = content
			if (@go_split_level2[0] eq " C")
			{
				$GO_all_C = $GO_all_C  . @go_split_level2[1] . "; ";
			}
			if (@go_split_level2[0] eq " P")
			{
				$GO_all_P = $GO_all_P  . @go_split_level2[1] . "; ";
			}
			if (@go_split_level2[0] eq " F")
			{
				$GO_all_F = $GO_all_F  . @go_split_level2[1] . "; ";
			}
		}
		if($_ =~m/^CC   -!- SUBCELLULAR LOCATION:/)
		{	
			@cc_split_level1=split(": ",$_);
			$output_loc =$output_loc . @cc_split_level1[1]; 
			$output_loc=~s/    / /g;
			$output_loc=~s/\W/ /g;
			$output_loc=~s/  / /g;
			$output_loc=~s/By similarity//g;	
		}
			if($_ =~m/^FT   SIGNAL/)
			{	
				$signal="Signal_peptide";
			}
	}	
	close tmp2;
	$quantile="n.d.";
	if ($Fc ne "n.d.")
	{		
		$quantile="middle";		
		if ($Fc < $low_percentile_cutoff)
		{
			$quantile="less";
		}
		if ($Fc > $high_percentile_cutoff)
		{
			$quantile="more";
		}
		if ($Fc < -3.5)
		{
			$h1++;
		}
		if (($Fc >= -3.5) && ($Fc < -2.5))
		{
			$h2++;
		}
		if (($Fc >= -2.5) && ($Fc < -1.5))
		{
			$h3++;
		}
		if (($Fc >= -1.5) && ($Fc < -0.5))
		{
			$h4++;
		}
		if (($Fc >= -0.5) && ($Fc < 0.5))
		{
			$h5++;
		}
		if (($Fc >= 0.5) && ($Fc < 1.5))
		{
			$h6++;
		}
		if (($Fc >= 1.5) && ($Fc < 2.5))
		{
			$h7++;
		}
		if (($Fc >= 2.5) && ($Fc < 3.5))
		{
			$h8++;
		}
		if ($Fc >= 3.5)
		{
			$h9++;
		}		
	}
	$joint_annotation = $GO_all_C. " " .$output_loc. " " .$signal;
	$sec = "no";
	#weak criteria first, these can be over-ruled by mitochondria
	if ((index($joint_annotation, "surface") != -1)|| (index($joint_annotation, "integral to membrane") != -1))
	{
		$sec = "yes";
	}
	if ((index($joint_annotation, "mitochondrion") != -1) || (index($joint_annotation, "mitochondrial") != -1)|| (index($joint_annotation, "nuclear membrane") != -1)|| (index($joint_annotation, "Nucleus membrane") != -1))
	{
		$sec = "no";
	}
	
	#now strong criteria which over-rule mitochondria
	if ((index($joint_annotation, "Signal_peptide") != -1) || (index($joint_annotation, "extracellular") != -1)|| (index($joint_annotation, "Secreted") != -1))
	{
		$sec = "yes";
	}
	$total_ids++;
	if($sec eq "yes")
	{
		$sec_ids++;
	}
	
	$output="yes";
	if(($filter_ntt eq "y" )&&($MS eq "")&&($DB eq ""))
	{
		$output="no";
	}
	$prime_tmp = $prime;
	$last_residue_of_prime = chop($prime_tmp);

	if($filter_ntt eq "y")
	{
		$output="no";
	}

	if(($filter_ntt eq "y" )&&($digestion_protease eq "t")&&(($last_residue_of_prime eq "R") ||($last_residue_of_prime eq "K")))
	{
		$output="yes";
	}
	if(($filter_ntt eq "y" )&&($digestion_protease eq "g")&&($last_residue_of_prime eq "E"))
	{
		$output="yes";
	}
	
	
	
	if($output eq "yes")
	{
		print outfile_non_red "$experiment\t$MS\t$DB\t$nonprime\t$prime\t$id1\t$id2\t$id3\t$description\t$Fc\t \t$quantile\t$Fc_orig\t$Position\t$Prot_Length\t$modifications\t$hyperscore\t$ppm\t$charge\t$mass\t$AltNames\t$GO_all_C\t$GO_all_P\t$GO_all_F\t$output_loc\t$signal\t$sec\n";
		$array_input=$Fc . "_".$Fc_error;
		push (@Fc_array,$array_input);	
	}	
}
open(merge_pdfs, ">merge_pdfs.bat");
print merge_pdfs "pdftk *.pdf cat output MS-MS_spectra.pdf\n del *light.pdf \n del *heavy.pdf\n del *acetylated.pdf\n ";
close merge_pdfs;
$command1= "mv merge_pdfs.bat ./pdf ";
$command2= "cd ./pdf && merge_pdfs.bat";
$command3="mv *TAILS_peptides_non-redundant.xls ./pdf";
close outfile_non_red;
system ($command1);
system ($command2);
system ($command3);
print "low cutoff = $low_percentile_cutoff   high cutoff = $high_percentile_cutoff\n";


#Fc histogram
my @data = (["<=-4", "-3", "-2", "-1", "0", "1", "2", "3", ">=4"],
            [$h1, $h2, $h3, $h4, $h5, $h6, $h7, $h8, $h9]);
my $mygraph = GD::Graph::bars->new(500, 500);
$mygraph->set(
    x_label     => 'Bins',
    y_label     => 'Number of IDs',
    title       => 'Fc-value distribution',
	x_label_position   => 0.5,
    transparent        => 0,
    bgclr              => 'white',
) or warn $mygraph->error;
my $myimage = $mygraph->plot(\@data) or die $mygraph->error;
$pic_file = $experiment."_histogramm.png";
$outfile = "histogramm.pdf";
open(PICTURE, ">$pic_file") or die("Cannot open file for writing");
# Make sure we are writing to a binary stream
binmode PICTURE;
# Convert the image to PNG and print it to the file PICTURE
print PICTURE $myimage->png;
close PICTURE;
system ("mv $pic_file ./pdf ");
#done Fc histogram

#localization pie chart
$other=$total_ids - $sec_ids;
my @data = (['sec. or surface', 'other'],
            [$sec_ids,$other]);
my $mygraph = GD::Graph::pie->new(500, 500);
$mygraph->set(
    title       => 'Secreted or cell surface associated',
    '3d'          => 1,
    transparent        => 0,
    bgclr              => 'white',
) or warn $mygraph->error;
$mygraph->set_value_font(GD::gdMediumBoldFont);
my $myimage = $mygraph->plot(\@data) or die $mygraph->error;
$pic_file = $experiment."_secreted.png";
open(PICTURE, ">$pic_file") or die("Cannot open file for writing");
# Make sure we are writing to a binary stream
binmode PICTURE;
# Convert the image to PNG and print it to the file PICTURE
print PICTURE $myimage->png;
close PICTURE;
system ("mv $pic_file ./pdf ");
#done localization pie chart


#localization hsitogram
my @data = (['sec. or surface', 'other'],
            [$sec_ids,$other]);
my $mygraph = GD::Graph::bars->new(500, 500);
$mygraph->set(
    x_label     => ' ',
    y_label     => 'Number of IDs',
    title       => 'Secreted or cell surface associated',
	show_values => 1,
	x_label_position   => 0.5,
    transparent        => 0,
    bgclr              => 'white',
) or warn $mygraph->error;
my $myimage = $mygraph->plot(\@data) or die $mygraph->error;
$pic_file = $experiment."_secreted-histogramm.png";
open(PICTURE, ">$pic_file") or die("Cannot open file for writing");
# Make sure we are writing to a binary stream
binmode PICTURE;
# Convert the image to PNG and print it to the file PICTURE
print PICTURE $myimage->png;
close PICTURE;
system ("mv $pic_file ./pdf ");
#done localization hsitogram
