## There are five lines of code that need to be modified before the example code can be used.  Each is marked with # <--

## Load the source code
source("C:\\Source code path\\Peaks_R_1.3.R")     # <--

## Load the C DLL ##
dll.path = "C:\\C_DLL_path\\Peaks_R121.dll"  	# <--
dyn.load(dll.path)  # Note: under Windows XP, using the compiled DLL provided, it is safe to ignore the warning about FPU control word


## Set the path of input files, the path of the metadata file, and the directory for output files ##
input.dir = "C:\\input_directory"		# <--
metadata.path = "C:\\metadata.txt"		# <--
output.dir = "C:\\output_directory"	# <--


## Call the peak extraction function. The default file.pattern matches all .csv files (.csv is the recommended format  ##
## for spectra exported from Ciphergen software). Parameters are described in comments of the getPeaks function.       ##
## The F test threshold used here by default is very strict, it should be adjusted to your individual need for quality ##
## vs number of peaks. Lower values of F.test.threshold will keep more peaks, lower values of chi.sq.threshold will    ##
## reject more spectra. ##

peaks <- getPeaks( spectra.dir=input.dir, files.pattern=".+\\.csv", metadata.path=metadata.path, additional.metadata=""  # give the input data
	, F.test.threshold=0.95, chi.sq.threshold=0.99   # set quality control thresholds
	, min.peak.widths=c(0.003, 0.008), max.peak.widths=c(0.01, 0.025), m.z.regions=c(3000, 18000)  # set peak width parameters
	, Memory.Conserve=FALSE, Peaks.Show=TRUE  # set the new flags introduced in version 1.2. The values here are the default values.
	, output.path.screened    = paste(output.dir, "peaks_screened.txt", sep="\\") # set output files ("" does not produce an output file)
	, output.path.all         = paste(output.dir, "peaks_all.txt"     , sep="\\")
	, output.path.averaged    = paste(output.dir, "peaks_averaged.txt", sep="\\")
	, output.path.peak.bounds = paste(output.dir, "peak_list.txt"     , sep="\\") )
	
	
## Perform Student's t test as an example of one way to identify differentially expressed peaks ##
	
if (peaks$return.code == 0) {
	print(noquote("Peak finding failed."))
} else {
	## T-test for differences between two groups, given by status1 and status2. The default values "WT" and "KO" are appropriate for the
	## example data.

	status1 = 3  # WT  R assigns numbers to each group in alphabetical order.
	status2 = 2  # KO  For data sets with only two groups, these values should be 1 and 2.

	p.threshold = 0.001  # Set the p value threshold for reporting differences.

	t.test.results = sapply( 1:dim(peaks$peaks.averaged)[[1]], FUN=function(index, data, status1, status2, status)
	{
		t.test(data[index, status==status1], data[index, status==status2])
	},peaks$peaks.averaged, status1, status2, peaks$status.averaged)

	for (index in (1:dim(t.test.results)[[2]])) {
		if (t.test.results["p.value", index] < p.threshold)
		{
			print( noquote( paste( "Peak ",  dimnames(peaks$peaks.averaged)[[1]][index]      , sep="" ) ) )
			print( noquote( paste( "T statistic = ", t.test.results["statistic", index][[1]] , sep="" ) ) )
			print( noquote( paste( "p value = ",     t.test.results["p.value",   index]      , sep="" ) ) )
			print( noquote( paste( status1, " mean = ",  t.test.results["estimate", index][[1]][1], sep="" ) ) )
			print( noquote( paste( status2, " mean = ",  t.test.results["estimate", index][[1]][2], sep="" ) ) )
			print( noquote( "" ) )
		}
	}
}