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


#MS2Oxoplot.pl sequentially reads in all .mgf-files in a user-defined directory/folder, extracts the intensity values of user-defined m/z values (+/- a user-defined mass-tolerance) and puts them out to .mgf-file specific .csv files (automated file name generation: INPUT_FILENAME_RepIntensities.csv).

#expected run-time on a "normal"desktop computer: < 60 sec.; increasing the number of reporter ions to extract (@unsortedRepIons) and the size of .mgf-file increases processing time.

#----User-definable PARAMETERS:

my @unsortedRepIons =(167.0914,183.0863,204.0867,222.0972,224.1118,243.0264,274.0921,284.0435,290.087,292.1027,308.0976,312.1289,316.1027,328.1238,332.0976,334.1133,350.1082,366.1395,370.1697,407.1661,495.1821,511.177,512.1974,528.1923,542.1716,553.224,569.21887,622.1284,657.2349,658.2553,673.2298,698.2615,699.2455,714.2564,715.2404,715.27677,792.3234,803.2928,819.2877,860.31427,876.30917,,948.3303,980.3201,446.09627,649.17567,731.27167,877.32957,1023.38747,1169.44537);

my $precision = 0.05; # +/- precision in amu

#----end of PARAMETERS

#----declaration of VARIABLES

my $MGF_OUTFILE_NAME;
my $REP_OUTFILE_NAME;
my $MGF_INFILE_NAME;
my $MASS;
my $TITLE;
my $PEPMASS; 
my $PRECURSOR_INTENSITY;
my $CHARGE;
my $RTINSECONDS = 0;
my $SCANS;
my $NUMBER_MS2SCANS_TOTAL;
my $NUMBER_MS2SCANS_FILTERED;

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

#----end of declaration of VARIABLES



my @RepIons = sort {$a <=> $b} @unsortedRepIons; #sort RepIons list by MoverZ

#----print PROMPT screen

print "---SugarQBits MS2Oxoplot.pl---\n";
print "User defined parameters:\n";
print "mass precision = +/-$precision amu\n";
print "Reporter Ions:\n";

foreach $x (@RepIons) {
	
	print $x.",";
}
print "\n---   ---\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);

#----endo of print PROMPT screen


#----creating OUTFILES

foreach $MGF_INFILE_NAME (@files) {

	$TITLE = $MGF_INFILE_NAME;
	
	$MGF_INFILE_NAME = join ("\\\\", $dirname, $MGF_INFILE_NAME);
	
	if ($MGF_INFILE_NAME =~ /.MGF/ || $MGF_INFILE_NAME =~ /.mgf/ ){
		
		print "SugarQbits RepX $MGF_INFILE_NAME...";

		open(INFILE,"<$MGF_INFILE_NAME") || die "$MGF_INFILE_NAME not found!\n";
		
		$NUMBER_MS2SCANS_TOTAL = 0; #reset counter for each new .mgf-file opened
		
		$NUMBER_MS2SCANS_FILTERED = 0;#reset counter for each new .mgf-file opened

		$REP_OUTFILE_NAME = ">>".substr ($MGF_INFILE_NAME,0,-4)."_RepIntensities.csv";

		open (OUTFILE2, $REP_OUTFILE_NAME);

		#---write header line for current .csv output-file

		print OUTFILE2 "MS2 Spectrum Title; Precursor mass; Charge;Retention Time";
				
		for ($x = 0; $x < @RepIons; $x++) { 
		
			print OUTFILE2 ";".$RepIons[$x];		
		
		}
		
		print OUTFILE2 ";BPI;TIC;Gscore\n";

#----end of creating OUTFILES


#---- READ IN DATA-FILE

		while (<INFILE>){

		chomp;
		

		@v = split (/=|\s/, $_);
	
		if ($v[0] eq "BEGIN"){
		
			#reset variables and empty arrays of data from previous MS2-spectrum in current .mgf-file
			
			@arr1 = ();
		
			@arr2 = ();
			
			@arr3 = ();
		
			$NUMBER_MS2SCANS_TOTAL++;
			
			$BPI = 0;
			
			$TIC = 0;
		
			#reset RepIons array arr2
			
			for ($x = 0; $x < @RepIons; $x++) { 
		
				$arr2[$x][0] = $RepIons[$x];
		
				$arr2[$x][1] = 0;
				
				$arr2[$x][2] = 0;
				
			
			}
									
		}elsif($v[0] eq "TITLE"){#MGF-MS2-spectrum title
				
			$TITLE = $v[1];	#DISCLAIMER: handles msconvert generated "TPP-compatible" spectrum titles; might truncate differntly constructed MS2-title formats
									
		}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{ #N.B. 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; #removes "+" or "-" sign from charge.
				
		}elsif ($v[0] eq "RTINSECONDS"){#MGF 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++;
		
			#--------DO ALL RepIntensities-file PROCESSING FROM HERE!				
			
			@arr3 = reverse sort { $a->[1] <=> $b->[1]} @arr1;
			
			
			for ($y = 0; $y < @arr3; $y++) { 
			
				for ($x = 0; $x < @RepIons; $x++) {
					
					if ($arr3[$y][0] < $arr2[$x][0] + $precision && $arr3[$y][0] > $arr2[$x][0] - $precision ){
						
						$arr2[$x][2] = $y+1;

					}

				}
					
			}
				
			print OUTFILE2 "TITLE=SQbRepX_".$TITLE.".".$NUMBER_MS2SCANS_FILTERED.";$PEPMASS;$CHARGE;$RTINSECONDS";
				
			for ($x = 0; $x < @RepIons; $x++) { 
		
				print OUTFILE2 ";".$arr2[$x][1];		
		
			}
				
			print OUTFILE2 ";".$BPI.";".$TIC."\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
		
			$v[0] = sprintf "%.6f",$v[0];
		
			$v[1] = sprintf "%.6f",$v[1];
		
			push @arr1, [$v[0], $v[1]];
			
			$TIC = $TIC + $v[1]; #TIC value is iteratively computed as sum of all intensity values in current MS2-spectrum
			
			if ($v[1] > $BPI) { #BPI value is iteratively picked as highest intensity value in current MS2-spectrum
			
				$BPI = $v[1];
			
			}
			
			if ($v[0] < $arr2[(@RepIons - 1)][0] + $precision){					
					
				for ($x = 0; $x < @RepIons; $x++) {
					
					if ($v[0] < $arr2[$x][0] + $precision && $v[0] > $arr2[$x][0] - $precision ){

						
						#find highest peak in mass-tolerance window
						
						if ($arr2[$x][1] < $v[1]) {

							$arr2[$x][1] = $v[1];

					
						}
					
					}

				}
					
			}
					
			
	
		}else{
		
			print "dont know what to do with: $v[0]-$v[1]\n"; #print out warning on unexpected .mgf-file format issues.
		
		}	
		
	}
	
	close OUTFILE2;

	close INFILE;
	
	print "done!\n($NUMBER_MS2SCANS_FILTERED RepXFiltered/$NUMBER_MS2SCANS_TOTAL total)\n";

}

}