#!/usr/local/bin/perl
use warnings;
use strict;
#use explicit;



#----PARAMETERS:

my $precision = 10; # precision in ppm
my $maxIsoptopeRatio = 10; # maximum difference in ratio of isoptope peak intensity;
my @unsortedGLYCAN_DB = ();
my $precursor_include = 0; #include propagation of original spectrum; "1" if to include;
my $number_ranks_reported = 1; #maximum number of ranks per MS2 spectrum propagated;
my $peptide_min_moverz = 600; #minimum moverz of peptide; parameter to avoid oxonium ions being wrongly selected as Y1.
my $trim_spectra = 1; #trims MS2 spectra after Y1 ion candidate (not the full MS2 spectrum goes to MGF file); anything else splits the MS2 data into peptide and glycan part.

#----end of PARAMETERS


my $MGF_OUTFILE_NAME;
my $Glyco_OUTFILE_NAME;
my $REP_OUTFILE_NAME;
my $MGF_INFILE_NAME;
my $MASS;
my $TITLE;
my $PEPMASS; 
my $PRECURSOR_INTENSITY;
my $CHARGE;
my $RTINSECONDS;
my $SCANS;
my $NUMBER_MS2SCANS_TOTAL;
my $NUMBER_MS2SCANS_FILTERED;
my $BPI;

my @arr1 = ();
my @arr2 = ();
my @v = ();
my $precursor;
my $x;
my $y;
my $i;	

print "---SugarQBits Kassonade---\n";
print "parameters:\n";
print "mass precision = +/-$precision ppm\n";
print "propagation of original MS2 spectrum (0=no): $precursor_include\n";
print "number of ranks propagated = $number_ranks_reported\n";
print "minimum peptide mass > $peptide_min_moverz amu\n";
print "trim MS2 spectra to [peptide+HexNAc]+ mass (1=yes): $trim_spectra\n";
print "---   ---\n";

print "INPUT (.mgf) FILE DIRECTORY (full path):>\n";

my $dirname = <>;

chop($dirname);

opendir(DIR, $dirname) or die "Could not open $dirname\n";

my @files = readdir(DIR);

shift @files; #get rid of "."

shift @files; #get rid of ".."

closedir(DIR);

foreach $MGF_INFILE_NAME (@files) {

	if ($MGF_INFILE_NAME =~ /.MGF/ || $MGF_INFILE_NAME =~ /.mgf/ ){
		
		print "SugarQbits Kassonade $MGF_INFILE_NAME...";

		$MGF_INFILE_NAME = join ("\\\\", $dirname, $MGF_INFILE_NAME);

		open(INFILE,"<$MGF_INFILE_NAME") || die "$MGF_INFILE_NAME not found!\n";

		$MGF_OUTFILE_NAME = ">>".substr ($MGF_INFILE_NAME,0,-4)."_Deglyco.mgf";
		
		open (OUTFILE, $MGF_OUTFILE_NAME);
		
		#$Glyco_OUTFILE_NAME = ">>".substr ($MGF_INFILE_NAME,0,-4)."_GlycoPeaks.mgf";
		
		#open (OUTFILE2, $Glyco_OUTFILE_NAME);

		#Start parsing MGF-file:

		while (<INFILE>){#feed in data line until EOF

			chomp;

			@v = split (/=|\s/, $_);
	
			if ($v[0] eq "BEGIN"){#MGF-MS2-spectrum begin
		
				#reset variables and empty arrays of data from previous MS2-spectrum
			
				@arr1 = ();
		
				@arr2 = ();
				
				$BPI = 0;
				
				$NUMBER_MS2SCANS_TOTAL++;
									
			}elsif($v[0] eq "TITLE"){#MGF-MS2-spectrum title
				
				$TITLE = $v[1];	#DISCLAIMER: will truncate msconvert generated "TPP-compatible" spectrum titles 
									
			}elsif ($v[0] eq "PEPMASS"){#MGF precursor ion m/z and its intensity
			
				$PEPMASS = $v[1];
			
				if (scalar @v > 2) {
			
					$PRECURSOR_INTENSITY = $v[2];
		
				}else{ #PEAKS generated mgf files do not contain precursor intensity values
		
					$PRECURSOR_INTENSITY = 0;
				}
				
			}elsif ($v[0] eq "CHARGE"){#MGF precursor charge state
				
				$CHARGE = substr $v[1], 0, 1; #remove + sign from charge.
				
			}elsif ($v[0] eq "RTINSECONDS"){#MGFclassics RetentionTimeINSECONDS
			
				$RTINSECONDS = $v[1];
			
			}elsif ($v[0] eq "SCANS"){#MGF MS2-spectrum scan number
			
				$SCANS = $v[1]
				
			}elsif ($v[0] eq "END"){
				
				$NUMBER_MS2SCANS_FILTERED++;
				
				#MGF file output
										
				$precursor = $PEPMASS * $CHARGE - ($CHARGE -1) * 1.00727;
					
				if ($precursor_include == 1) {
					
					unshift @arr2,[$precursor,$BPI]; #adds the original spectrum precursor intensity of BasePeak. relevant Y1 peaks should be higher? 
				
				} 

				@v = reverse sort {$a->[1] <=> $b->[1]}@arr2; # sort @arr2 by descending intensity
				
				@arr2 = @v; #copy @v back into @arr2

				my $i = 0;
				
				for ($x = 0; $x < @arr2; $x++) {
				
					if ($arr2[$x][1] > 0 && $i < $number_ranks_reported){
					
						$i++;
						
						if (($precursor - $arr2[$x][0]) == 0) {
							
							print OUTFILE "BEGIN IONS\nTITLE=SQbitK_".$TITLE.":".$i.":Gdelta~0"."\nPEPMASS=".($arr2[$x][0])." $PRECURSOR_INTENSITY\nCHARGE=1+\n";
							
						}else{
						
							print OUTFILE "BEGIN IONS\nTITLE=SQbitK_".$TITLE.":".$i.":Gdelta~".($precursor - $arr2[$x][0] + 203.0794)."\nPEPMASS=".($arr2[$x][0])." $PRECURSOR_INTENSITY\nCHARGE=1+\n";

						}

						#print OUTFILE2 "BEGIN IONS\n";
					
						#print OUTFILE2 "TITLE=SQbitK_".$TITLE.":".$i.":Gdelta~".($precursor - $arr2[$x][0] + 203.0794); #print complete glycan delta mass
						
						#print OUTFILE2 "\n";

						#print OUTFILE2 "PEPMASS=".($arr2[$x][0])." $PRECURSOR_INTENSITY\n"; #set potential Y1 as new precursor mass 
							
						#print OUTFILE2 "CHARGE=1+\n";
						
						for ($y = 0 ; $y < @arr1; $y++){
				
							if ($trim_spectra eq 1 && $arr1[$y][1] > 0 && $arr1[$y][0] < ($arr2[$x][0] + 203.0794 + 1 )) { #only print items with intensity > 0 and lower then the new precursor
					
								print OUTFILE $arr1[$y][0]." ".$arr1[$y][1]."\n";
					
							#}elsif ($trim_spectra == 1 && $arr1[$y][1] > 0 && $arr1[$y][0] > ($arr2[$x][0] + 203.0794 + 1) ) { #only print items with intensity > 0 and lower then the new precursor
					
								#print OUTFILE2 ($arr1[$y][0]-$arr2[$x][0])." ".$arr1[$y][1]."\n";
					
							}elsif ($trim_spectra ne 1 && $arr1[$y][1] > 0){
						
								print OUTFILE $arr1[$y][0]." ".$arr1[$y][1]."\n";
					
							}
				
						}
					
						print OUTFILE "END IONS\n";
						
						#print OUTFILE2 "END IONS\n";
						
					}
				
				}
						
			}elsif ($v[0] !~ /[a-df-zA-DF-Z]/ && $v[1] !~ /[a-df-zA-DF-Z]/) { #Lines that do not contain letters (but "E" for exponential numbers) are deemed MS2 data pairs
				
				push @arr1,[$v[0], $v[1]]; #push and store current values into @arr1
				
				if ($BPI < $v[1]) { $BPI = $v[1]; } #find BasePeak of current MS2 scan;
				
				#check for glycodistances to previously imported moverz values
				
				if ($v[0] > $peptide_min_moverz	 && $v[1] > 0){
				
					for ($x = @arr1 - 2; $x > 0; $x --){

						if (($v[0] - $arr1[$x][0]) < 120.042259 + ($v[0]/1000000)*$precision && ($v[0] - $arr1[$x][0]) > 120.042259 - ($v[0]/1000000)*$precision) {								
								
							$v[1] = $v[1] * $arr1[$x][1]; #reward N-glyco specific loss for Y1-ion picking by inflating intensity value put into @arr2
															
						}elsif (($v[0] - $arr1[$x][0]) < 146.0579 + ($v[0]/1000000)*$precision && ($v[0] - $arr1[$x][0]) > 146.0579 - ($v[0]/1000000)*$precision) {								
								
							$v[1] = $v[1] / $arr1[$x][1]; #punish dHex specific loss for Y1-ion picking by decreasing intensity value put into @arr2
															
						}elsif (($v[0] - $arr1[$x][0]) < 162.0528 + ($v[0]/1000000)*$precision && ($v[0] - $arr1[$x][0]) > 162.0528 - ($v[0]/1000000)*$precision) {								
								
							$v[1] = $v[1] / $arr1[$x][1]; #punish Hex specific loss for Y1-ion picking by decreasing intensity value put into @arr2
															
						}elsif (($v[0] - $arr1[$x][0]) < 203.0794 + ($v[0]/1000000)*$precision && ($v[0] - $arr1[$x][0]) > 203.0794 - ($v[0]/1000000)*$precision) {								
								
							push @arr2, [$v[0],$v[1]]; #the current moverz is put onto potential Y1-ion list for picking: ranking is based on pseudo-Intensity
															
						}elsif (($v[0] - $arr1[$x][0]) > 203.0794 + ($v[0]/1000000)*$precision){
						
							$x = 0; #stop looping
						}
					
					}
			
				}
			
			}else{
			
				print "dont know what to do with: $v[0]-$v[1]\n";

			}			
		
		}
	
	close OUTFILE;
	
	#close OUTFILE2;

	close INFILE;

	print "done!\n($NUMBER_MS2SCANS_FILTERED GFiltered/$NUMBER_MS2SCANS_TOTAL total)\n";

	}

}