## ----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()