Mercurial > repos > pieterlukasse > prims_metabolomics

## read args:
args <- commandArgs(TRUE)
# xset data:
args.xsetData <- args[1]

args.rtStart  <- strtoi(args[2])
args.rtEnd <- strtoi(args[3])

args.mzStart <- as.double(args[4])
args.mzEnd <- as.double(args[5])
# there are 2 options: specify a mz range or a mz list:
if (args.mzStart < 0)
{
	args.mzList <- as.double(strsplit(args[6], ",")[[1]])
	cat(typeof(as.double(strsplit(args[6], ",")[[1]])))
	args.mzTolPpm <- as.double(args[7])
	# calculate mzends based on ppm tol:
	mzListEnd <- c()
	mzListStart <- c()
	for (i in 1:length(args.mzList))
	{
		mzEnd <- args.mzList[i] + args.mzList[i]*args.mzTolPpm/1000000.0
		mzStart <- args.mzList[i] - args.mzList[i]*args.mzTolPpm/1000000.0
		mzListEnd <- c(mzListEnd, mzEnd)
		mzListStart <- c(mzListStart, mzStart)
	}
	str(mzListStart)
	str(mzListEnd)
} else {
	mzListEnd <- c(args.mzEnd)
	mzListStart <- c(args.mzStart)
}

args.sampleNames <- strsplit(args[8], ",")[[1]]
# trim leading and trailing spaces:
args.sampleNames <- gsub("^\\s+|\\s+$", "", args.sampleNames)

args.combineSamples <- args[9]
args.rtPlotMode <- args[10]

## report files
args.htmlReportFile <- args[11]
args.htmlReportFile.files_path <- args[12]


if (length(args) == 13)
{
	args.outLogFile <- args[13]
	# suppress messages:
	# Send all STDERR to STDOUT using sink() see http://mazamascience.com/WorkingWithData/?p=888
	msg <- file(args.outLogFile, open="wt")
	sink(msg, type="message")
	sink(msg, type="output")
}

# TODO - add option to do masses in same plot (if given in same line oid) or in separate plots
# TODO2 - let it run in parallel

tryCatch(
        {
	        library(metaMS)

			# load the xset data :
			xsetData <- readRDS(args.xsetData)
			# if here to support both scenarios:
			if ("xcmsSet" %in% slotNames(xsetData) )
			{
				xsetData <- xsetData@xcmsSet
			}

			# report
			dir.create(file.path(args.htmlReportFile.files_path), showWarnings = FALSE, recursive = TRUE)
			message(paste("\nGenerating report.........in ", args.htmlReportFile.files_path))

			html <- "<html><body><h1>Extracted Ion Chromatograms (EIC) matching criteria</h1>"

			if (args.combineSamples == "No")
			{
				if (length(args.sampleNames) > 1 && length(mzListStart) > 1 && length(args.sampleNames) != length(mzListStart))
					stop(paste("The number of sample names should match the number of m/z values in the list. Found ", length(mzListStart),
					          " masses while ",  length(args.sampleNames), " sample names were given."))

		  		iterSize <- length(args.sampleNames)
				# these can be set to 1 or 0 just as a trick to iterate OR not over the items. If the respective list is of length 1, only the first item should be used
				fixSampleIdx <- 1
				fixMzListIdx <- 1
				if (length(args.sampleNames) == 1)
				{
					fixSampleIdx <- 0
					iterSize <- length(mzListStart)
				}
				if (length(mzListStart) == 1)
				{
					fixMzListIdx <- 0
				}
				lineColors <- rainbow(iterSize)
				for (i in 0:(iterSize-1))
				{
					message("\nGetting EIC... ")
					eiccor <- getEIC(xsetData,
										mzrange=matrix(c(mzListStart[i*fixMzListIdx+1],mzListEnd[i*fixMzListIdx+1]),nrow=1,ncol=2,byrow=TRUE),
										rtrange=matrix(c(args.rtStart,args.rtEnd),nrow=1,ncol=2,byrow=TRUE),
										sampleidx=c(args.sampleNames[i*fixSampleIdx+1]), rt=args.rtPlotMode)

					message("\nPlotting figures... ")
					figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="")
					html <- paste(html,"<img src='", "figure", i,".png' /><br/>", sep="")
					png( figureName, type="cairo", width=1100,height=250 )
					#plot(eiccor, col=lineColors[i+1])
					# black is better in this case:
					plot(eiccor)
					legend('topright', # places a legend at the appropriate place
							legend=c(args.sampleNames[i*fixSampleIdx+1]), # puts text in the legend
							lty=c(1,1), # gives the legend appropriate symbols (lines)
							lwd=c(2.5,2.5))

					devname = dev.off()
				}

			} else {
				for (i in 1:length(mzListStart))
				{
					message("\nGetting EIC... ")
					eiccor <- getEIC(xsetData,
										mzrange=matrix(c(mzListStart[i],mzListEnd[i]),nrow=1,ncol=2,byrow=TRUE),
										rtrange=matrix(c(args.rtStart,args.rtEnd),nrow=1,ncol=2,byrow=TRUE),
										sampleidx=args.sampleNames, rt = args.rtPlotMode)

										#set size, set option (plot per sample, plot per mass)

					message("\nPlotting figures... ")
					figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="")
					html <- paste(html,"<img src='", "figure", i,".png' />", sep="")
					png( figureName, type="cairo", width=1100,height=450 )
					lineColors <- rainbow(length(args.sampleNames))
					plot(eiccor, col=lineColors)
					legend('topright', # places a legend at the appropriate place
						  legend=args.sampleNames, # puts text in the legend
						  lty=c(1,1), # gives the legend appropriate symbols (lines)
						  lwd=c(2.5,2.5),
				          col=lineColors)
					devname = dev.off()
				}
			}
			if (args.rtPlotMode == "corrected")
			{
				html <- paste(html,"<p>*rt values are corrected ones</p></body><html>")
			}
			html <- paste(html,"</body><html>")
			message("finished generating report")
			write(html,file=args.htmlReportFile)
			# unlink(args.htmlReportFile)
			cat("\nWarnings================:\n")
			str( warnings() )
		},
        error=function(cond) {
            sink(NULL, type="message") # default setting
			sink(stderr(), type="output")
            message("\nERROR: ===========\n")
            print(cond)
        }
    )
author	pieter.lukasse@wur.nl
date	Thu, 11 Dec 2014 14:32:19 +0100
parents	f772a5caa86a
children