6
+ − 1 # shows all alignment results in a rt region
+ − 2 ## read args:
+ − 3 args <- commandArgs(TRUE)
+ − 4 # xset data:
+ − 5 args.xsetData <- args[1]
+ − 6
+ − 7 args.rtStart <- strtoi(args[2])
+ − 8 args.rtEnd <- strtoi(args[3])
+ − 9
+ − 10 # limit max diff to 600 and minNrSamples to at least 2 :
+ − 11 if (args.rtEnd - args.rtStart > 600)
+ − 12 stop("The RT window should be <= 600")
+ − 13
+ − 14 args.minNrSamples <- strtoi(args[4]) #min nr samples
+ − 15 if (args.minNrSamples < 2)
+ − 16 stop("Set 'Minimum number of samples' to 2 or higher")
+ − 17
+ − 18
+ − 19 args.sampleNames <- strsplit(args[5], ",")[[1]]
+ − 20 # trim leading and trailing spaces:
+ − 21 args.sampleNames <- gsub("^\\s+|\\s+$", "", args.sampleNames)
+ − 22
+ − 23 ## report files
+ − 24 args.htmlReportFile <- args[6]
+ − 25 args.htmlReportFile.files_path <- args[7]
+ − 26
+ − 27
+ − 28 if (length(args) == 8)
+ − 29 {
+ − 30 args.outLogFile <- args[8]
+ − 31 # suppress messages:
+ − 32 # Send all STDERR to STDOUT using sink() see http://mazamascience.com/WorkingWithData/?p=888
+ − 33 msg <- file(args.outLogFile, open="wt")
+ − 34 sink(msg, type="message")
+ − 35 sink(msg, type="output")
+ − 36 }
+ − 37
+ − 38
+ − 39
+ − 40 tryCatch(
+ − 41 {
+ − 42 library(metaMS)
+ − 43
+ − 44 # load the xset data :
+ − 45 xsetData <- readRDS(args.xsetData)
+ − 46 # if here to support both scenarios:
+ − 47 if ("xcmsSet" %in% slotNames(xsetData) )
+ − 48 {
+ − 49 xsetData <- xsetData@xcmsSet
+ − 50 }
+ − 51
+ − 52 # report
+ − 53 dir.create(file.path(args.htmlReportFile.files_path), showWarnings = FALSE, recursive = TRUE)
+ − 54 message(paste("\nGenerating report.........in ", args.htmlReportFile.files_path))
+ − 55
+ − 56 write(paste("<html><body><h1>Extracted Ion Chromatograms (EIC) of alignments with peaks in ",args.minNrSamples, " or more samples</h1>"),
+ − 57 file=args.htmlReportFile)
+ − 58
+ − 59 gt <- groups(xsetData)
+ − 60 message("\nParsed groups... ")
+ − 61 groupidx1 <- which(gt[,"rtmed"] > args.rtStart & gt[,"rtmed"] < args.rtEnd & gt[,"npeaks"] >= args.minNrSamples) # this should be only on samples selected
+ − 62
+ − 63 if (length(groupidx1) > 0)
+ − 64 {
+ − 65 message("\nGetting EIC... ")
+ − 66 eiccor <- getEIC(xsetData, groupidx = c(groupidx1), sampleidx=args.sampleNames)
+ − 67 #eicraw <- getEIC(xsetData, groupidx = c(groupidx1), sampleidx=args.sampleNames, rt = "raw")
+ − 68
+ − 69 #sampleNamesIdx <- which(sampnames(LC$xset@xcmsSet) %in% args.sampleNames, arr.ind = TRUE)
+ − 70 #or (from bioconductor code for getEIC): NB: this is assuming sample indexes used in data are ordered after the order of sample names!!
+ − 71 sampNames <- sampnames(xsetData)
+ − 72 sampleNamesIdx <- match( args.sampleNames, sampNames)
+ − 73 message(paste("Samples: ", sampleNamesIdx))
+ − 74
+ − 75 #TODO html <- paste(html, "<table><tbody>")
+ − 76 message(paste("\nPlotting figures... "))
+ − 77
+ − 78 #get the mz list (interestingly, this [,"mz"] is relatively slow):
+ − 79 peakMzList <- xsetData@peaks[,"mz"]
+ − 80 peakSampleList <- xsetData@peaks[,"sample"]
+ − 81 #signal to noise list:
+ − 82 peakSnList <- xsetData@peaks[,"sn"]
+ − 83
+ − 84 message(paste("Total nr of peaks: ", length(peakMzList)))
+ − 85
+ − 86 for (i in 1:length(groupidx1))
+ − 87 {
+ − 88 groupMembers <- xsetData@groupidx[[groupidx1[i]]]
+ − 89
+ − 90 groupMzList <- ""
+ − 91 groupSampleList <- ""
+ − 92 finalGroupSize <- 0
+ − 93
+ − 94 for (j in 1:length(groupMembers))
+ − 95 {
+ − 96 #get only the peaks from the selected samples:
+ − 97 memberSample <- peakSampleList[groupMembers[j]]
+ − 98 memberSn <- peakSnList[groupMembers[j]]
+ − 99 if (!is.na(memberSn) && memberSample %in% sampleNamesIdx)
+ − 100 {
+ − 101 message(paste("Group: ", groupidx1[i], " / Member sample: ", memberSample))
+ − 102 memberMz <- peakMzList[groupMembers[j]]
+ − 103
+ − 104
+ − 105 if (finalGroupSize == 0)
+ − 106 {
+ − 107 groupMzList <- memberMz
+ − 108 groupSampleList <- sampNames[memberSample]
+ − 109 } else {
+ − 110 groupMzList <- paste(groupMzList,",",memberMz, sep="")
+ − 111 groupSampleList <- paste(groupSampleList,",",sampNames[memberSample], sep="")
+ − 112 }
+ − 113
+ − 114 finalGroupSize <- finalGroupSize +1
+ − 115 }
+ − 116 }
+ − 117 # here we do the real check on group size and only the groups that have enough
+ − 118 # members with signal to noise > 0 will be plotted here:
+ − 119 if (finalGroupSize >= args.minNrSamples)
+ − 120 {
+ − 121 message(paste("Plotting figure ",i, " of max ", length(groupidx1)," figures... "))
+ − 122
+ − 123 figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="")
+ − 124 write(paste("<img src='", "figure", i,".png' />", sep="") ,
+ − 125 file=args.htmlReportFile, append=TRUE)
+ − 126
+ − 127 png( figureName, type="cairo", width=800 )
+ − 128 plot(eiccor, xsetData, groupidx = i)
+ − 129 devname = dev.off()
+ − 130
+ − 131 write(paste("<p>Alignment id: [", groupidx1[i], "]. m/z list of peaks in alignment with signal/noise <> NA:<br/>", groupMzList,"</p>", sep="") ,
+ − 132 file=args.htmlReportFile, append=TRUE)
+ − 133 write(paste("<p>m/z values found in the following samples respectively: <br/>", groupSampleList,"</p>", sep="") ,
+ − 134 file=args.htmlReportFile, append=TRUE)
+ − 135 }
+ − 136 }
+ − 137
+ − 138 }
+ − 139
+ − 140 write("</body><html>",
+ − 141 file=args.htmlReportFile, append=TRUE)
+ − 142 message("finished generating report")
+ − 143 # unlink(args.htmlReportFile)
+ − 144 cat("\nWarnings================:\n")
+ − 145 str( warnings() )
+ − 146 },
+ − 147 error=function(cond) {
+ − 148 sink(NULL, type="message") # default setting
+ − 149 sink(stderr(), type="output")
+ − 150 message("\nERROR: ===========\n")
+ − 151 print(cond)
+ − 152 }
+ − 153 )