diff --git a/patRoon_codes/Suspect_pos_neg.R b/patRoon_codes/Suspect_pos_neg.R
new file mode 100644
index 0000000000000000000000000000000000000000..021f03ad0f01df747a5cb41d0fe95f91517e8bf7
--- /dev/null
+++ b/patRoon_codes/Suspect_pos_neg.R
@@ -0,0 +1,111 @@
+################################Suspect screening of CSF with patRoon 2.1.0###############################
+
+###Load library
+library(patRoon)
+
+
+###Set all options
+#Set the path of the working directory
+setwd("../../../../patRoon_codes/")
+
+#set the name of the output folder
+output <- "report_suspect"
+
+# Modify this value to point to the location of the MetFrag2.4.5-CL.jar file
+options(patRoon.path.MetFragCL = "~/opt/metfrag/MetFrag2.4.5-CL.jar")
+
+options(patRoon.path.obabel = "~/babel/bin/")
+
+options(patRoon.path.MetFragPubChemLite = ".../patRoon_files/databases/PubChemLite_exposomics_20220729.csv")
+
+
+#Load the csv with the data of the mzML files (note that here only the RP files are shown,as an example but there are two different analysis.csv files por the HILIC analysis)
+#Positive ionization files
+anaInfoPos <- read.csv("analysis_RP_POS.csv")
+
+#Negative ionization files
+anaInfoNeg <- read.csv("analysis_RP_NEG.csv")
+
+
+#Select the path of the suspect list
+suspFile <- read.csv(".../patRoon_files/suspects/name_of_the_suspect_list.csv")
+
+
+###Extract features
+library(xcms)
+
+#To find positive features
+fListPos <- findFeatures(anaInfoPos, "xcms3", param = xcms::CentWaveParam(ppm = 5, peakwidth = c(10, 30), mzdiff = -0.001, mzCenterFun = "wMean", integrate = 1L, snthresh = 10, prefilter = 
+                                                                            c(3, 100), noise = 7500, fitgauss = FALSE, firstBaselineCheck = TRUE),verbose = TRUE)
+#To find negative features
+
+fListNeg <- findFeatures(anaInfoNeg, "xcms3", param = xcms::CentWaveParam(ppm = 5, peakwidth = c(10, 30), mzdiff = -0.001, mzCenterFun = "wMean", integrate = 1L, snthresh = 10, prefilter = 
+                                                                            c(3, 100), noise = 7500, fitgauss = FALSE, firstBaselineCheck = TRUE),verbose = TRUE)
+#To merge pos and neg results
+fList <- makeSet(fListPos, fListNeg, adducts = c("[M+H]+", "[M-H]-"))
+
+
+###Group features
+library(BiocParallel)
+BiocParallel::register(BiocParallel::SerialParam(), default = TRUE)
+
+fGroups <- groupFeatures(fList, "xcms3")
+
+library ("MetaClean") 
+fGroups <- filter(fGroups, absMinIntensity = 100,
+                  relMinReplicateAbundance = NULL, maxReplicateIntRSD = NULL,
+                  blankThreshold = 3, removeBlanks = TRUE,
+                  retentionRange = NULL, mzRange = NULL)
+
+###Screen suspects
+fGroups <- screenSuspects(fGroups, suspFile, rtWindow = 12, mzWindow = 0.005, onlyHits = TRUE)
+
+###Retrieve MS peak lists
+avgPListParams <- getDefAvgPListParams(clusterMzWindow = 0.005)
+mslists <- generateMSPeakLists(fGroups, "mzr", maxMSRtWindow = 5, 
+                               precursorMzWindow = 0.5, avgFeatParams = 
+                               avgPListParams, avgFGroupParams =
+                               avgPListParams)
+
+mslists <- filter(mslists, withMSMS = TRUE, absMSIntThr = NULL, absMSMSIntThr = NULL, relMSIntThr = NULL,
+                        relMSMSIntThr = NULL, topMSPeaks = NULL, topMSMSPeaks = 25,
+							deIsotopeMS = FALSE, deIsotopeMSMS = FALSE)
+	
+							 
+###Find compound structure candidates 
+compounds <- generateCompounds(fGroups, mslists, "metfrag", method = "CL",
+                               dbRelMzDev = 5, fragRelMzDev = 5, 
+                               fragAbsMzDev = 0.002, 
+                               database = "pubchemlite",
+                               scoreTypes = c("fragScore","metFusionScore","score", "individualMoNAScore",
+                                              "AnnoTypeCount","Patent_Count", "DisorderDisease", "PubMed_Count"),
+                               maxCandidatesToStop = 100, timeoutRetries = 20)  
+							   
+###Annotate suspects
+fGroups <- annotateSuspects(fGroups, formulas = NULL,
+								compounds = compounds, MSPeakLists = mslists)
+
+###Report results
+
+#Create the csv tables with groups and peak intensities for the set (positive and negative results combined)
+reportCSV(fGroups, path = output, reportFeatures = FALSE,
+          compounds = compounds, compoundsNormalizeScores = "max")
+		  
+#Individual features tables 
+pos_features <- as.data.table(unset(fGroups, "positive")) # table with all feature group data for positive mode
+write.csv(pos_features, "report_suspect/pos_features.csv")
+neg_features <- as.data.table(unset(fGroups, "negative")) #  table with all feature group data for negative mode
+write.csv(neg_features, "report_suspect/neg_features.csv")
+
+
+# Summary of MetFrag Results in a a Single Table 
+MFsummary <- as.data.table(compounds)
+outputSummary <- paste(output, "MFsummary.csv", sep = "/")
+write.csv(MFsummary, outputSummary)
+
+#Plots of every group of features (in the set)
+reportPDF(fGroups, path = output, reportFGroups = TRUE, 
+          reportFormulaSpectra = TRUE, compounds = compounds, 
+          compoundsNormalizeScores = "max",
+          MSPeakLists = mslists)
+