## ----LibraryPreload, message=FALSE--------------------------------------- library(Risa) library(xcms) library(CAMERA) library(pcaMethods) ## ----settings------------------------------------------------------------ ## How many CPU cores has your machine (or cluster) ? nSlaves=1 # prefilter <- c(3,200) ## standard prefilter=c(6,750) ## quick-run for debugging ## ----rISA, cache=TRUE---------------------------------------------------- ISAmtbls2 <- readISAtab(find.package("mtbls2")) a.filename <- ISAmtbls2["assay.filenames"][[1]] ## ----PeakPicking, cache=TRUE, warning=FALSE------------------------------ mtbls2Set <- processAssayXcmsSet(ISAmtbls2, a.filename, method="centWave", prefilter=prefilter, snthr=25, ppm=25, peakwidth=c(5,12), nSlaves=nSlaves) ## ----xcmsSet------------------------------------------------------------- show(mtbls2Set) ## ----retcor-------------------------------------------------------------- mtbls2Set <- group(mtbls2Set, minfrac=1, bw=3) retcor(mtbls2Set, plottype="mdevden") ## ----QCintensity--------------------------------------------------------- boxplot(groupval(mtbls2Set, value="into")+1, col=as.numeric(sampclass(mtbls2Set))+1, log="y", las=2) ## ----fillPeaks----------------------------------------------------------- mtbls2Set <- fillPeaks(mtbls2Set, nSlaves=nSlaves) ## ----plotQC-------------------------------------------------------------- plotQC(mtbls2Set) ## ----QCPCA, fig.show='hold'---------------------------------------------- sdThresh <- 4.0 ## Filter low-standard deviation rows for plot data <- log(groupval(mtbls2Set, value="into")+1) pca.result <- pca(data, nPcs=3) plotPcs(pca.result, type="loadings", col=as.numeric(sampclass(mtbls2Set))+1) ## ----CAMERA, warning=FALSE, results='hide'------------------------------- an <- xsAnnotate(mtbls2Set, sample=seq(1,length(sampnames(mtbls2Set))), nSlaves=nSlaves) an <- groupFWHM(an) an <- findIsotopes(an) # optional but recommended. an <- groupCorr(an, graphMethod="lpc", calcIso = TRUE, calcCiS = TRUE, calcCaS = TRUE, cor_eic_th=0.5) ## Setup ruleSet rs <- new("ruleSet") rs@ionlistfile <- file.path(find.package("mtbls2"), "lists","ions.csv") rs@neutraladditionfile <- file.path(find.package("mtbls2"), "lists","neutraladdition.csv") rs@neutrallossfile <- file.path(find.package("mtbls2"), "lists","neutralloss.csv") rs <- readLists(rs) rs <- setDefaultParams(rs) rs <- generateRules(rs) an <- findAdducts(an, rules=rs@rules, polarity="positive") ## ----diffreport---------------------------------------------------------- dr <- diffreport(mtbls2Set, sortpval=FALSE, filebase="mtbls2diffreport", eicmax=20 ) cspl <- getPeaklist(an) annotatedDiffreport <- cbind(dr, cspl) ## ----diffreportPspec----------------------------------------------------- interestingPspec <- tapply(seq(1, nrow(annotatedDiffreport)), INDEX=annotatedDiffreport[,"pcgroup"], FUN=function(x, a) {m <- median(annotatedDiffreport[x, "pvalue"]); p <- max(annotatedDiffreport[x, "pcgroup"]); as.numeric(c(pvalue=m,pcgroup=p))}, annotatedDiffreport) interestingPspec <- do.call(rbind, interestingPspec) colnames(interestingPspec) <- c("pvalue", "pcgroup") o <- order(interestingPspec[,"pvalue"]) pdf("interestingPspec.pdf") dummy <- lapply(interestingPspec[o[1:40], "pcgroup"], function(x) {suppressWarnings(plotPsSpectrum(an, pspec=x, maxlabel=5))}) dev.off() ## ----assembleMAF--------------------------------------------------------- pl <- annotatedDiffreport charge <- sapply(an@isotopes, function(x) { ifelse( length(x) > 0, x$charge, NA) }) abundance <- groupval(an@xcmsSet, value="into") ## ## load ISA assay files ## a.samples <- ISAmtbls2["samples.per.assay.filename"][[ a.filename ]] ## ## These columns are defined by mzTab ## maf.std.colnames <- c("identifier", "chemical_formula", "description", "mass_to_charge", "fragmentation", "charge", "retention_time", "taxid", "species", "database", "database_version", "reliability", "uri", "search_engine", "search_engine_score", "modifications", "smallmolecule_abundance_sub", "smallmolecule_abundance_stdev_sub", "smallmolecule_abundance_std_error_sub") ## ## Plus the columns for the sample intensities ## all.colnames <- c(maf.std.colnames, a.samples) ## ## Now assemble new maf ## l <- nrow(pl) maf <- data.frame(identifier = character(l), chemical_formula = character(l), description = character(l), mass_to_charge = pl$mz, fragmentation = character(l), charge = charge, retention_time = pl$rt, taxid = character(l), species = character(l), database = character(l), database_version = character(l), reliability = character(l), uri = character(l), search_engine = character(l), search_engine_score = character(l), modifications = character(l), smallmolecule_abundance_sub = character(l), smallmolecule_abundance_stdev_sub = character(l), smallmolecule_abundance_std_error_sub = character(l), abundance, stringsAsFactors=FALSE) ## ----exportMAF----------------------------------------------------------- ## ## Make sure maf table is quoted properly, ## and add to the ISAmtbls2 assay file. ## maf_character <- apply(maf, 2, as.character) write.table(maf_character, file=paste(tempdir(), "/a_mtbl2_metabolite profiling_mass spectrometry_maf.csv", sep=""), row.names=FALSE, col.names=all.colnames, quote=TRUE, sep="\t", na="\"\"") ISAmtbls2 <- updateAssayMetadata(ISAmtbls2, a.filename, "Metabolite Assignment File", paste(tempdir(), "/a_mtbl2_metabolite profiling_mass spectrometry_maf.csv", sep="")) write.assay.file(ISAmtbls2, a.filename) ## ----sessionInfo--------------------------------------------------------- sessionInfo()