DIMaSP
Direct Infusion Mass Spectrometry Processing

Step 2: Main processing
Release Version 1.0

Copyright © 2017 Arthur T Zielinski, Stephen J Fuller, Claudio Bortolini, Ivan Kourtchev, Markus Kalberer

This file is part of “Direct Infusion Mass Spectrometry Processing (DIMaSP)”.
Direct Infusion Mass Spectrometry Processing (DIMaSP) is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. 
Direct Infusion Mass Spectrometry Processing (DIMaSP) is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with Direct Infusion Mass Spectrometry Processing (DIMaSP).  If not, see <http://www.gnu.org/licenses/>.

This notebook will take the paired sample and blank files listed in tablename located in the specified directory, folder, and process them based on (the following is a partial list):
	1. Remove entries below a specified noise (intensity) level unique to each sample and blank
	2. Filter peaks in the sample that have a corresponding match in the blank with an intensity above the Minimum Sample-to-Blank Ratio (specified in tablename)
	3. Remove entries from the “clean” sample file with drifts (i.e. Delta (ppm)) outside of the range set by the Sample upper ppm	limit and Sample lower ppm limit in tablename.
	4. Remove entries without a Carbon or Hydrogen
	5. The processed Hydrogen count depends on polarity:
		ESI NEG: H+1
		ESI POS: H if contains Na, otherwise (H-1)
		APPI POS: H if integer RDB, otherwise (H-1)
		APPI NEG: H if integer RDB, otherwise (H+1)
	6. Check H/C, N/C, O/C, and S/C ratios based on the specified limits below (e.g. hcmin). Note: H/C ratio depends on the source and polarity (e.g. ESI Negative)
	7. Remove entries with non-integer double bound equivalent (DBE)
	8. Remove entries that fail the Nitrogen Rule
	9. Remove entries that contain C13/S34 if there is an equivalent peak without them. This process accounts for the isotopic ratios for Carbon and Sulphur (see “Other Options” below). The isotopic removal process does not function for assignments with both C13 and S34.

The script can process files that used ESI Negative, ESI Positive, APPI Negative, or APPI Positive as specified in tablename.

Each sample and blank combination will be output as a separate .csv file named based on the sample name with the processed mode appended (i.e. “-esi-neg”, “-esi-pos”, “-appi-neg”, or “-appi-pos”). If duplicates have been deleted, the file will also have “-deldup” appended to the name. The output file will the original data as well the adjusted data (based on source/polarity). Additional columns include:
	1. Abs. Corrected Delta: the absolute value of Delta (ppm) after being corrected for the mass shift. Used to filter duplicates.
	2. O/C and H/C: the value of the modified O/C and H/C ratios, respectively.
	3. OSC, KM, and KMD: the Carbon Oxidation State, Kendrick Mass, and Kendrick Mass Defect, respectively.

To run the code go to Evaluation > Evaluate Notebook. Code will be completed once “Running...” no longer appears at the top of the window. While the code is running it will display the entry in tablename it is processing. After each file it will output the number of assignments (including duplicate peaks) after various stages of processing. Upon completion it will display the time taken to process all of the files in tablename.

Example Input:
	folder = “C:\\Users\\Name\\Raw”;
	tablename = “MSP_Input.csv”;
	
Other Options:
	deleteduplicates is a toggle (1 == yes, 0 == no) for whether or not to delete duplicate peaks. The peak with the minimum 	Abs. Corrected Delta is kept.
	hcmin, hcmax, ocmax, ocmin, ncmax, scmax are the ratio limits (min and max) for H/C, O/C, N/C, and S/C
	upper13c, lower13c, upper34s, lower34s are toggles (0 or 1) for checking isotopic ratios. They should remain as having 	the “lower” values set to 0, and the “upper” values set to 1.
	c1213 isotopic ratio (for check) of 12C/13C. Default 0.011.
	s3234 isotopic ratio (for check) of 32S/34S. Default 0.045.
	scalec13 isotopic ratio scaling factor of 12C/13C. Default 1.0.
	scales34 isotopic ratio scaling factor of 12C/13C. Default 1.0.

Important notes:
	1. The main body of the code is hidden to prevent any accidental changes but can be accessed by selecting (left click) the 	hidden cell (short line on right of notebook) and going to Cell > Cell Properties > Open and Cell > Cell Properties > Editable.

Inputs

folder="C:\\Users\\Name\\Raw"; 
(* Folder containing the samples, blanks, and table file *)
tablename="MassSpecProcessingTable.csv"; 
(* Table filename INCLUDING file extension (i.e. ".csv"); The table must contain : the sample name, sample mass shift, Sample SD, Number of compound used to SDcalculation for the sample, Samble mass drift filter factor,sample noise, blank name, Blank mass shift, Blank SD, Number of compound used to SDcalculation for the Blank, Blank noise, Sample to blank ratio, acquisition mode*)
deleteduplicates=0; (* Toggle for deleting duplicates. 1 = Yes, 0 = No *)
hcmin=0.3; (* H/C minimum *)
hcmax=2.5; (* H/C maximum *)
ocmax=2; (* O/C maximum *)
ocmin=0; (* O/C minimum *)
ncmax=1.3; (* N/C maximum *)
scmax=0.8; (* S/C maximum *)
upper13c=1; (* DO NOT ADJUST *)
lower13c=0; (* DO NOT ADJUST *)
upper34s=1; (* DO NOT ADJUST *)
lower34s=0; (* DO NOT ADJUST *)
c1213 = 0.011; (* 12C/13C Ratio *)
s3234=0.045; (* 32S/34S Ratio *)
scalec13=1.0; (* Scaling factor for 12C/13C ratio *)
scales34=1.0; (* Scaling factor for 32S/34S ratio *)

Reading information from table

SetDirectory[folder]; (* Change directory to folder *)
speclist=Drop[Import[tablename],1]; (* Table containing specifications for each file imported from table with headers removed *)
filename=Table[speclist[[i,1]] <> ".csv",{i,Dimensions[speclist][[1]]}]; (* Column of sample filenames with .csv added *)
sampledrift=speclist[[All,2]]; (* Column of sample mass drifts *)
samplestdev=speclist[[All,3]]; (* Column of sample standard deviation *)
stdev=Sqrt[((speclist[[All,3]])^2)+((speclist[[All,10]])^2)]; (* Column of root mean square of sample and blank standard deviations *)
df=speclist[[All,4]]+speclist[[All,11]]-2; (* Column of degrees of freedoms (for two sample test) *)
tstud=Table[0,{i,Length[speclist]}]; (* Initialise tstud column *)
Do[
tstud[[i]]= If[df[[i]]>30,1.96,Switch[df[[i]],1,12.71,2,4.303,3,3.182,4,2.776,5,2.571,6,2.447,7,2.365,8,2.306,9,2.262,10,2.228,11, 2.201,12,2.179, 13,2.160,14,
2.145,15,2.131,16,2.120,17,2.110,18,2.101,19,2.093,20,2.086,21,2.080,22,2.074,23,2.069,24,2.064,25,2.060,
26,2.056,27,2.052,28,2.048,29,2.045,30,2.042]],
{i,Length[speclist]}]; (* Student's t-Distribution for confidence level of 95% *)
(*sampledriftfactor=speclist[[All,5]]; (* Column of multiplicative factors for drift filtering (e.g. 2) *) (* CAN REMOVE *)*)
ppmupper=speclist[[All,6]]; (* Column of sample ppm upper limit for mass drift filter *)
ppmlower=speclist[[All,7]]; (* Column of sample ppm lower limit for mass drift filter *)
samplenoise=speclist[[All,5]]; (* Column of sample noise limits *)
blankname=Table[speclist[[i,8]] <> ".csv",{i,Dimensions[speclist][[1]]}]; (* Column of blank filenames with .csv added *)
blankdrift=speclist[[All,9]]; (* Column of blank mass drifts *)
blanknoise=speclist[[All,12]]; (* Column of blank noise limits *)
snb=speclist[[All,13]]; (* Column of sample to blank intensity ratios *)
sourcepolarity=speclist[[All,14]]; (* Column of source (APPI/ESI) and polarity (POS/NEG) distinction *)

Blank Subtraction and Filtering

TimeSt=AbsoluteTime[]; (* Start Timer *)
Do[ (* Outer loop goes through each sample and blank combination *)
(* Print details of current run *)
Print["Entry " <> ToString[f]];If[Min[df[[f]]]≤0,Print["Error: Mass Drift does not have enough degrees of freedom"]];
Print["Mode: " <> sourcepolarity[[f]]];
Print["Sample: " <> filename[[f]]];
Print["Sample Drift: " <> ToString[sampledrift[[f]]]];
Print["Blank: " <> blankname[[f]]];
Print["Blank Drift: " <> ToString[blankdrift[[f]]]];
Print["t-student: " <> ToString[tstud[[f]]]];
Print["Combined SD: " <> ToString[stdev[[f]]]<> "\n"];

rawsamplewithheaders=Import[folder<>"\\"<>filename[[f]]]; (* Sample data with headers *)
rawblankwithheaders=Import[folder<>"\\"<>blankname[[f]]]; (* Blank data with headers *)
rawsample=rawsamplewithheaders[[9;;Length[rawsamplewithheaders]]];  (* Sample data *)
stg0=rawsample; (* For checking # of rows later *)
rawblank=rawblankwithheaders[[9;;Length[rawblankwithheaders]]]; (* Blank data *)
Do[ (* Duplicate masses and intensities as needed - SAMPLE *)
		prow=rawsample[[i-1]]; (* Previous Row in Sample data *)
		row=rawsample[[i]]; (* Current Row in Sample data *)
		If[NumberQ[row[[1]]]==False,rawsample[[i,1;;3]]=prow[[1;;3]]] (* If first entry in current row is NOT a number, then set the first three entries of current row to be the same as previous row *)
		,{i,2,Length[rawsample]} (* Loop through entire sample *)
];
Do[ (* Duplicate masses and intensities as needed - BLANK *)
		prow=rawblank[[i-1]]; (* Same as above, but for blank *)
		row=rawblank[[i]];
		If[NumberQ[row[[1]]]==False,rawblank[[i,1;;3]]=prow[[1;;3]]]
		,{i,2,Length[rawblank]}
];
sample1 = Select[rawsample,(#[[2]]>samplenoise[[f]])&&(NumberQ[#[[4]]])&]; (* Remove sample data with intensities less than noise, and any empty compounds *)
blank1=Select[rawblank,(#[[2]]>blanknoise[[f]])&&(NumberQ[#[[4]]])&]; (* Remove blank data with intensities less than noise, and any empty compounds *)
(* NEW, v1.37 *)
stg1=sample1;  (* For checking # of rows later *)

newcols={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; (* zeros for new columns *)
samplec13s34check=Table[{sample1[[i]],newcols},{i,Length[sample1]}]; (* Add extra 0 columns to data *)
samplec13s34check=Partition[Flatten[samplec13s34check],31]; (* Changed all partitions to be 30 from 25 for new cols *)
Do[
form=StringSplit[samplec13s34check[[i,7]]];
Do[
elem=Characters[form[[j]]];
If[Length[elem]>1,If[elem[[2]]=="a",elem={"Y"}]];
If[Length[elem]>4,If[elem[[5]]=="C",elem={"X"}]];
If[Length[elem]>4,If[elem[[5]]=="S",elem={"Z"}]];
elemL=Length[elem];
elemN=If[elemL==1,1,If[elemL==2,ToExpression[elem[[2]]],If[elemL==3,(ToExpression[elem[[2]]]*10+ToExpression[elem[[3]]]),(ToExpression[elem[[2]]]*100+ToExpression[elem[[3]]]*10+ToExpression[elem[[4]]])]]]; (* Changed to allow for 100s of element *)
If[elem[[1]]=="N",samplec13s34check[[i,10]]=elemN,
If[elem[[1]]=="C",samplec13s34check[[i,8]]=elemN,
If[elem[[1]]=="H",samplec13s34check[[i,9]]=elemN,
If[elem[[1]]=="O",samplec13s34check[[i,12]]=elemN,
If[elem[[1]]=="S",samplec13s34check[[i,11]]=elemN,
If[elem[[1]]=="Y",samplec13s34check[[i,15]]=elemN,
If[elem[[1]]=="X",samplec13s34check[[i,13]]=elemN,
If[elem[[1]]=="Z",samplec13s34check[[i,14]]=elemN
]]]]]]]];
,{j,Length[form]}];
,{i,Length[samplec13s34check]}];
(* /NEW, v1.37 *)

sample2=sample1; (* Copy sample data *)
sample2[[All,1]]=sample2[[All,1]]-sampledrift[[f]]*sample2[[All,1]]/1000000; (* Add drift *)
blank2=blank1; (* Copy blank data *)
blank2[[All,1]]=blank2[[All,1]]-blankdrift[[f]]*blank2[[All,1]]/1000000; (* Add drift *)
sample2cols=DeleteDuplicates[sample2[[All,1;;2]]]; (* Extract first two columns and delete duplicates - NOTE: only deletes consecutive duplicates (I think) *)
blank2cols=DeleteDuplicates[blank2[[All,1;;2]]]; (* Extract first two columsn and delete duplicates - NOTE: only deletes consecutive duplicates (I think) *)
samplecheck=Table[0,{Length[sample2cols]}]; (* Empty array for use in next loop *)
Do[
blankcheck=Table[3,{Length[blank2cols]}];
Do[
(* If |Ms-Mb|/(|Ms+Mb|/2)*1000000 < 2*1.96 stdev AND Is/Ib > Sample to Blank Intensity Ratio, then add a counter to "samplecheck" array in corresponding row *) (* Replaced 1.96 with tstud *)
If[(Abs[sample2cols[[scounter,1]]-blank2cols[[bcounter,1]]]/(
Abs[sample2cols[[scounter,1]]+blank2cols[[bcounter,1]]]/2)*1000000>2 * tstud[[f]] * stdev[[f]]), (* Check if blank is outside standard deviation of sample entry *)
blankcheck[[bcounter]]=1, (* If outside, set check to 1 *)
If[sample2cols[[scounter,2]]/blank2cols[[bcounter,2]]>snb[[f]], (* If not, check if the sample-to-blank ratio is above the specified limit *)
blankcheck[[bcounter]]=1, (* If yes, set check to 1 *)
blankcheck[[bcounter]]=0]] (* If no, set check to 0 --> To be deleted *)
,{bcounter,1,Length[blank2cols]}
];
samplecheck[[scounter]]=Count[blankcheck,0] (* Reverse count (want to count instances when there are NO matches from above loop *)
, {scounter,1,Length[sample2cols]}
];
samplecolscomb=Partition[Flatten[Transpose[{sample2cols,samplecheck}]],3]; (* Combine samplecheck with sample2cols*)
goodsamplemass=Select[samplecolscomb,#[[3]]==0&]; (* Only select rows in combined array with non-zero samplecheck entries (i.e. entries that passed above test) *)
sampleBS=Select[sample2,Count[goodsamplemass[[All,1]],#[[1]]]>0&];  (* From "big" sample table, only select those rows that have a matching mass in the goodsamplemass array - this is the blank substracted table *)
stg2=sampleBS;  (* For checking # of rows later *)
xcal=sampleBS; (* Continue on from blank substracted data *)
(* New in v1.3 *)
xcal[[All,1]]=Round[xcal[[All,1]]+sampledrift[[f]]*xcal[[All,1]]/1000000,0.00001]; (* Remove drift, Round to same digits as before *)

(*titles=Partition[Flatten[{rawsamplewithheaders[[7]],"C","H","N","S","O","13C","34S","Na","Theo N.Mass","DBE","N.Comp","C","H","N","S","O","13C","34S", "O/C", "H/C", "OSC", "KM", "KMD"}],30];*)
titles=Flatten[{rawsamplewithheaders[[7]],"C","H","N","S","O","13C","34S","Na","Theo N.Mass","DBE","N.Comp","C","H","N","S","O","13C","34S", "Abs. Corrected Delta","O/C", "H/C", "OSC", "KM", "KMD"}];

Do[
prow=xcal[[i-1]];
row=xcal[[i]];
If[NumberQ[row[[1]]]==False,xcal[[i,1;;3]]=prow[[1;;3]]];,{i,2,Length[xcal]}];
(*=Select[xcal,NumberQ[#[[4]]]==True&];*) (* Not needed since already removed above *)
xcal1=xcal;
newcols={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; (* zeros for new columns *)
xcal2=Table[{xcal1[[i]],newcols},{i,Length[xcal1]}]; (* Add extra 0 columns to data *)
xcal3=Partition[Flatten[xcal2],31]; (* Changed all partitions to be 30 from 25 for new cols *)

Do[
form=StringSplit[xcal3[[i,7]]];
Do[
elem=Characters[form[[j]]];
If[Length[elem]>1,If[elem[[2]]=="a",elem={"Y"}]];
If[Length[elem]>4,If[elem[[5]]=="C",elem={"X"}]];
If[Length[elem]>4,If[elem[[5]]=="S",elem={"Z"}]];
elemL=Length[elem];
elemN=If[elemL==1,1,If[elemL==2,ToExpression[elem[[2]]],If[elemL==3,(ToExpression[elem[[2]]]*10+ToExpression[elem[[3]]]),(ToExpression[elem[[2]]]*100+ToExpression[elem[[3]]]*10+ToExpression[elem[[4]]])]]]; (* Changed to allow for 100s of element *)
If[elem[[1]]=="N",xcal3[[i,10]]=elemN,
If[elem[[1]]=="C",xcal3[[i,8]]=elemN,
If[elem[[1]]=="H",xcal3[[i,9]]=elemN,
If[elem[[1]]=="O",xcal3[[i,12]]=elemN,
If[elem[[1]]=="S",xcal3[[i,11]]=elemN,
If[elem[[1]]=="Y",xcal3[[i,15]]=elemN,
If[elem[[1]]=="X",xcal3[[i,13]]=elemN,
If[elem[[1]]=="Z",xcal3[[i,14]]=elemN
]]]]]]]];
,{j,Length[form]}];
,{i,Length[xcal3]}];

xcal4a=xcal3;
(*(* Claudio method for checking drift *)xcal4=Select[xcal4a,#[[5]]≥(sampledrift[[f]]-sampledriftfactor[[f]] samplestdev[[f]])&&#[[5]]<=(sampledrift[[f]]+sampledriftfactor[[f]] samplestdev[[f]])&]; (* Newly changed... *)*)
xcal4=Select[xcal4a,#[[5]]≥ppmlower[[f]]&&#[[5]]≤ppmupper[[f]]&];
stg3=xcal4;  (* For checking # of rows later *)
xcal5a=Select[xcal4,(#[[8]]+#[[13]])>0&];
xcal5=Select[xcal5a,#[[9]]>0&];
(* Added switch statement for polarity *)
xcal6a=Switch[sourcepolarity[[f]], (* Select H/C check based on mode *)
"ESI NEG",Select[xcal5,(#[[9]]+1)/(#[[8]]+#[[13]])<=hcmax&],
"ESI POS", Select[xcal5,(#[[9]]-If[#[[15]]==0,1,0])/(#[[8]]+#[[13]])<=hcmax&],
"APPI POS", Select[xcal5,(#[[9]]-If[IntegerQ[#[[6]]],0,1])/(#[[8]]+#[[13]])<=hcmax&],
"APPI NEG", Select[xcal5,(#[[9]]+If[IntegerQ[#[[6]]],0,1])/(#[[8]]+#[[13]])<=hcmax&]];
xcal6b=Switch[sourcepolarity[[f]], (* Select H/C check based on mode *)
"ESI NEG",Select[xcal6a,(#[[9]]+1)/(#[[8]]+#[[13]])>=hcmin&],
"ESI POS", Select[xcal6a,(#[[9]]-If[#[[15]]==0,1,0])/(#[[8]]+#[[13]])>=hcmin&],
"APPI POS", Select[xcal6a,(#[[9]]-If[IntegerQ[#[[6]]],0,1])/(#[[8]]+#[[13]])>=hcmin&],
"APPI NEG", Select[xcal6a,(#[[9]]+If[IntegerQ[#[[6]]],0,1])/(#[[8]]+#[[13]])>=hcmin&]];
(*xcal6a=Select[xcal5,(#[[9]]-If[#[[15]]<1,1,0])/(#[[8]]+#[[13]])<=hcmax&];
xcal6b=Select[xcal6a,(#[[9]]-If[#[[15]]<1,1,0])/(#[[8]]+#[[13]])>=hcmin&];*)
xcal6c=Select[xcal6b,#[[10]]/(#[[8]]+#[[13]])<=ncmax&];
xcal6d=Select[xcal6c,#[[12]]/(#[[8]]+#[[13]])<=ocmax&];
xcal6d2=Select[xcal6d,#[[12]]/(#[[8]]+#[[13]])>=ocmin&];
xcal6=Select[xcal6d2,#[[11]]/(#[[8]]+#[[13]])<=scmax&];
stg4=xcal6;  (* For checking # of rows later *)
Do[
xcal6[[k,19]]=xcal6[[k,8]];
xcal6[[k,21;;25]]=xcal6[[k,10;;14]];
xcal6[[k,20]]=Switch[sourcepolarity[[f]], (* Select processing method for Hydrogen column *)
"ESI NEG",xcal6[[k,9]]+1,
"ESI POS", If[xcal6[[k,15]]==0,xcal6[[k,9]]-1,xcal6[[k,9]]],
"APPI POS", If[IntegerQ[xcal6[[k,6]]],xcal6[[k,9]],xcal6[[k,9]]-1],
"APPI NEG", If[IntegerQ[xcal6[[k,6]]],xcal6[[k,9]],xcal6[[k,9]]+1]];
(*If[IntegerQ[xcal6[[k,6]]],xcal6[[k,20]]=xcal6[[k,9]],xcal6[[k,20]]=xcal6[[k,9]]-1];*)
xcal6[[k,18]]="C"<>ToString[xcal6[[k,19]]]<>
If[xcal6[[k,24]]<1,"",If[xcal6[[k,24]]<2,"[13]C","[13]C"<>ToString[xcal6[[k,24]]]]]<>
"H"<>ToString[xcal6[[k,20]]]<>
If[xcal6[[k,21]]<1,"",If[xcal6[[k,21]]<2,"N","N"<>ToString[xcal6[[k,21]]]]]<>
If[xcal6[[k,22]]<1,"",If[xcal6[[k,22]]<2,"S","S"<>ToString[xcal6[[k,22]]]]]<>
If[xcal6[[k,25]]<1,"",If[xcal6[[k,25]]<2,"[34]S","[34]S"<>ToString[xcal6[[k,25]]]]]<>
If[xcal6[[k,23]]<1,"",If[xcal6[[k,23]]<2,"O","O"<>ToString[xcal6[[k,23]]]]];
xcal6[[k,17]]=xcal6[[k,19]]-xcal6[[k,20]]/2+xcal6[[k,21]]/2+1;
xcal6[[k,16]]=((xcal6[[k,19]]*12)+(xcal6[[k,20]]*1.007825)+(xcal6[[k,21]]*14.003074)+(xcal6[[k,22]]*31.972072)+(xcal6[[k,23]]*15.994915)+(xcal6[[k,24]]*13.003355)+(xcal6[[k,25]]*33.967868));
xcal6[[k,26]]=Abs[xcal6[[k,5]]-sampledrift[[f]]]; (* Adjust mass drift for sample ABSOLUTE VALUE *)
xcal6[[k,27]]=N[xcal6[[k,23]]/xcal6[[k,19]]]; (* O/C Ratio *)
xcal6[[k,28]]=N[xcal6[[k,20]]/xcal6[[k,19]]]; (* H/C Ratio *)
xcal6[[k,29]]=2 xcal6[[k,27]]-xcal6[[k,28]]; (* O_SC *)
xcal6[[k,30]]=xcal6[[k,16]] (14/14.01565); (* KM *)
xcal6[[k,31]]=N[12 xcal6[[k,19]]+xcal6[[k,20]]+14 xcal6[[k,21]]+16 xcal6[[k,23]]+32 xcal6[[k,22]]-xcal6[[k,30]]]; (* KMD *)
,{k,Length[xcal6]}];

xcal7=Select[xcal6,(#[[17]]-Round[#[[17]]])==0&&#[[17]]>=0&]; (* DBE filtering -> Removes any non-integer DBE (and any negative DBE) *)
output=Table[{},{tt,Length[xcal7]+1}];
xcal8a=Sort[xcal7,#1[[2]]>#2[[2]]&];
xcal8=Select[xcal8a,EvenQ[(((#[[19]]+#[[24]])*12)+#[[20]]+(#[[21]]*14)+((#[[22]]+#[[25]])*32)+(#[[23]]*16)+#[[21]])]==True&];
xcal8=SortBy[xcal8,{-#[[2]],#[[26]]}&]; (* Re-sort data based on adjusted intensity AND mass delta *)
stg5=xcal8;  (* For checking # of rows later *)

output[[1]]=titles;
(*null={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
c12s32only=Select[xcal8,(#[[24]]+#[[25]])==0&];
c13ors34=Select[xcal8,(#[[24]]+#[[25]])==1&];
c13only=Select[c13ors34,#[[24]]==1&];
s34only=Select[c13ors34,#[[25]]==1&];
c13checka:=Select[c12s32only,(#[[19]]-c13only[[kk,19]])==1&];
c13checkb:=Select[c13checka,(#[[20]]-c13only[[kk,20]])==0&];
c13checkc:=Select[c13checkb,(#[[23]]-c13only[[kk,23]])==0&];
c13checkd:=Select[c13checkc,(#[[21]]-c13only[[kk,21]])==0&];
c13checke:=Select[c13checkd,(#[[22]]-c13only[[kk,22]])==0&];
c13checkf:=Select[c13checke,(c13only[[kk,24]]-#[[24]])==1&];
c13checkI:=Select[c13checkf,(c13only[[kk,2]]/#[[2]])<(c1213*#[[8]]*upper13c)&&(c13only[[kk,2]]/#[[2]])>(c1213*#[[8]]*lower13c)&];
c13founda=Table[{If[Length[c13checkI]==0,null,c13only[[kk]]]},{kk,Length[c13only]}];
c13foundb=Partition[Flatten[c13founda],31];
c13found=Select[c13foundb,#[[1]]>0&];
s34checka:=Select[c12s32only,(#[[19]]-s34only[[kk,19]])==0&];
s34checkb:=Select[s34checka,(#[[20]]-s34only[[kk,20]])==0&];
s34checkc:=Select[s34checkb,(#[[21]]-s34only[[kk,21]])==0&];
s34checkd:=Select[s34checkc,(#[[22]]-s34only[[kk,22]])==1&];
s34checke:=Select[s34checkd,(#[[23]]-s34only[[kk,23]])==0&];
s34checkf:=Select[s34checke,(s34only[[kk,25]]-#[[25]])==1&];
s34checkI:=Select[s34checkf,(s34only[[kk,2]]/#[[2]])<(s3234*#[[11]]*upper34s)&&(s34only[[kk,2]]/#[[2]])>(s3234*#[[11]]*lower34s)&];
s34founda=Table[{If[Length[s34checkI]==0,null,s34only[[kk]]]},{kk,Length[s34only]}];
s34foundb=Partition[Flatten[s34founda],31];
s34found=Select[s34foundb,#[[1]]>0&];
searchc13:=Select[c13found,(#[[1]]-c12s32only[[gg,1]])==0&];
c13removeda=Table[{If[Length[searchc13]==0,c12s32only[[gg]],null]},{gg,Length[c12s32only]}];
c13removedb=Partition[Flatten[c13removeda],31];
c13removed=Select[c13removedb,#[[1]]>0&];
searchs34:=Select[s34found,(#[[1]]-c13removed[[gg,1]])==0&];
s34removeda=Table[{If[Length[searchs34]==0,c13removed[[gg]],null]},{gg,Length[c13removed]}];
s34removedb=Partition[Flatten[s34removeda],31];
s34removed=Select[s34removedb,#[[1]]>0&];
c13s34removed=Select[s34removed,(#[[13]]+#[[14]])==0&];*)
null={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
c12s32onlybig=Select[samplec13s34check,(#[[13]]+#[[14]])==0&];
c12s32only=Select[xcal8,(#[[13]]+#[[14]])==0&];
c13ors34=Select[xcal8,(#[[13]]+#[[14]])==1&];
c13only=Select[c13ors34,#[[13]]==1&];
s34only=Select[c13ors34,#[[14]]==1&];
c13checka:=Select[c12s32onlybig,(#[[8]]-c13only[[kk,8]])==1&];
c13checkb:=Select[c13checka,(#[[9]]-c13only[[kk,9]])==0&];
c13checkc:=Select[c13checkb,(#[[12]]-c13only[[kk,12]])==0&];
c13checkd:=Select[c13checkc,(#[[10]]-c13only[[kk,10]])==0&];
c13checke:=Select[c13checkd,(#[[11]]-c13only[[kk,11]])==0&];
c13checkf:=Select[c13checke,(c13only[[kk,13]]-#[[13]])==1&];
(* Pre v1.40 
c13checkI:=Select[c13checkf,(c13only[[kk,2]]/#[[2]])<(c1213*#[[8]]*upper13c)&&(c13only[[kk,2]]/#[[2]])>(c1213*#[[8]]*lower13c)&]; *)
(* New in v1.40 *)
c13checkI:=Select[c13checkf,(c13only[[kk,2]]/#[[2]])<(scalec13*c1213*#[[8]]*upper13c)&&(c13only[[kk,2]]/#[[2]])>(scalec13*c1213*#[[8]]*lower13c)&];
(* /New *)
c13founda=Table[{If[Length[c13checkI]==0,null,c13only[[kk]]]},{kk,Length[c13only]}];
c13foundb=Partition[Flatten[c13founda],31];
c13found=Select[c13foundb,#[[1]]>0&];
s34checka:=Select[c12s32onlybig,(#[[8]]-s34only[[kk,8]])==0&];
s34checkb:=Select[s34checka,(#[[9]]-s34only[[kk,9]])==0&];
s34checkc:=Select[s34checkb,(#[[10]]-s34only[[kk,10]])==0&];
s34checkd:=Select[s34checkc,(#[[11]]-s34only[[kk,11]])==1&];
s34checke:=Select[s34checkd,(#[[12]]-s34only[[kk,12]])==0&];
s34checkf:=Select[s34checke,(s34only[[kk,14]]-#[[14]])==1&];
s34checkI:=Select[s34checkf,(s34only[[kk,2]]/#[[2]])<(scales34*s3234*#[[11]]*upper34s)&&(s34only[[kk,2]]/#[[2]])>(scales34*s3234*#[[11]]*lower34s)&];
s34founda=Table[{If[Length[s34checkI]==0,null,s34only[[kk]]]},{kk,Length[s34only]}];
s34foundb=Partition[Flatten[s34founda],31];
s34found=Select[s34foundb,#[[1]]>0&];
(* Change again *)
searchc13:=Select[c13found,(#[[1]]-c12s32only[[gg,1]])==0&];
c13removeda=Table[{If[Length[searchc13]==0,c12s32only[[gg]],null]},{gg,Length[c12s32only]}];
c13removedb=Partition[Flatten[c13removeda],31];
c13removed=Select[c13removedb,#[[1]]>0&];
searchs34:=Select[s34found,(#[[1]]-c13removed[[gg,1]])==0&];
s34removeda=Table[{If[Length[searchs34]==0,c13removed[[gg]],null]},{gg,Length[c13removed]}];
s34removedb=Partition[Flatten[s34removeda],31];
s34removed=Select[s34removedb,#[[1]]>0&];
c13s34removed=Select[s34removed,(#[[13]]+#[[14]])==0&];



output2=Table[{},{tt,Length[c13s34removed]+1}];
output2[[1]]=titles;

(*Do[
output[[m+1]]=xcal8[[m]];,{m,Length[xcal8]}]; (* Unuesd output *)*)

Do[
output2[[m+1]]=c13s34removed[[m]];
,{m,Length[c13s34removed]}];



filenameappend=Switch[sourcepolarity[[f]], (* Select file naming system *)
"ESI NEG","-esi-neg.csv",
"ESI POS", "-esi-pos.csv",
"APPI POS", "-appi-pos.csv",
"APPI NEG", "-appi-neg.csv"];

(* Export data. Delete duplicates if desired *)
If[deleteduplicates==0,Export[StringDrop[filename[[f]],-4]<>filenameappend,output2 ],If[deleteduplicates==1, Export[StringDrop[filename[[f]],-4]<>"-deldup"<>filenameappend,DeleteDuplicates[output2,#1[[1]]==#2[[1]]&] ],Print["Error: Invalid delete duplicates value. Must be 0 or 1."]]];

(* DEBUGGING *)
d0=Dimensions[stg0]; (* RAW *)
d1=Dimensions[stg1]; (* Removing blank formulae and noise *)
d2=Dimensions[stg2]; (* Blank subtracted *)
d3=Dimensions[stg3]; (* ppm upper/ ppm lower *)
d4=Dimensions[stg4]; (* Carbon ratio checks *)
d5=Dimensions[stg5]; (* DBE & Nitrogen rule *)
d6=Dimensions[output2]-1; (* Final (isotopes). Minus one because it counts header *)
Print["Raw: ",d0[[1]],"\nNoise (and empties) removed: ",d1[[1]],"\nBlank subtracted: ",d2[[1]], "\nppm upper/lower filtering: ",d3[[1]], "\nCarbon ratio checks: ",d4[[1]], "\nDBE and Nitrogen Rule: ",d5[[1]], "\nFinal (isotopes): ",d6[[1]]]


(* \DEBUGGING *)

,{f,Length[filename]}]

TimeEn=AbsoluteTime[];
Print["Total Processing Time (minutes): " <> ToString[(TimeEn-TimeSt)/60]];