Faked xcms diffreport

fakediffreport.R

library(multtest)


## A convenience function cut&paste from
## https://github.com/sneumann/xcms/blob/master/R/xcmsSet.R#L2052
pval <- function(X, classlabel, teststat) {

    n1 <- rowSums(!is.na(X[,classlabel == 0]))
    n2 <- rowSums(!is.na(X[,classlabel == 1]))
    A <- apply(X[,classlabel == 0], 1, sd, na.rm=TRUE)^2/n1 ## sd(t(X[,classlabel == 0]), na.rm = TRUE)^2/n1
    B <- apply(X[,classlabel == 1], 1, sd, na.rm=TRUE)^2/n2 ## sd(t(X[,classlabel == 1]), na.rm = TRUE)^2/n2
    df <- (A+B)^2/(A^2/(n1-1)+B^2/(n2-1))

    pvalue <- 2 * (1 - pt(abs(teststat), df))
    invisible(pvalue)
}


## Fake some data
dr <- data.frame(mz=100:120,
                 rt=10:30,
                 tstat=NA,
                 pvalue=NA,
                 fold=NA,
                 anova=NA,                 
                 samp1a=rnorm(21, mean=17),
                 samp1b=rnorm(21, mean=17),
                 samp1c=rnorm(21, mean=17),
                 samp1d=rnorm(21, mean=17),
                 samp2a=rnorm(21, mean=38),
                 samp2b=rnorm(21, mean=38),
                 samp2c=rnorm(21, mean=38),
                 samp2d=rnorm(21, mean=38),
                 samp3a=rnorm(21, mean=24),
                 samp3b=rnorm(21, mean=24),
                 samp3c=rnorm(21, mean=24),
                 samp3d=rnorm(21, mean=24))

## T-Test between class samp1 and samp2

testval <- dr[, grep("samp1|samp2", colnames(dr))]
testclab <- as.factor(substr(colnames(testval), 1, 5))

dr[,"tstat"] <- mt.teststat(testval, testclab) ## well, actually: , ...)
dr[,"pvalue"] <- pval(testval, as.integer(as.factor(testclab))-1, tstat)

## Plain foldchange between class samp1 and samp2
dr[,"fold"] <- ( rowMeans(dr[, grep("samp1", colnames(dr))])
    / rowMeans(dr[, grep("samp2", colnames(dr))]) )

## Do the anova stuff:
values <- dr[, grep("samp", colnames(dr))]

pvalAnova<-c()
for(i in 1:nrow(dr)){
    var <- as.numeric(values[i,])
    ano <- summary(aov(var ~ substr(colnames(values), 1, 5)))
    pvalAnova <- append(pvalAnova, unlist(ano)["Pr(>F)1"])
}

dr[,"anova"] <- pvalAnova

## Now annotate the columns:
attr(dr$pvalue, "stato") <- "[stato;4711;t-test]"