MS Differential Expression Tests

### R code from vignette source 'msmsTests-Vignette2.Rnw'

###################################################
### code chunk number 1: msmsTests-Vignette2.Rnw:49-50
###################################################
options(continue=" ")


###################################################
### code chunk number 2: SetupChunk
###################################################
library(msmsTests)
data(msms.dataset)
msms.dataset
msms.counts <- exprs(msms.dataset)
dim(msms.counts)
table(pData(msms.dataset)$treat,pData(msms.dataset)$batch)


###################################################
### code chunk number 3: DistChunk
###################################################
idx <- grep("HUMAN",featureNames(msms.dataset)) 
mSpC <- t( apply(msms.counts[idx,],1,function(x) 
                   tapply(x,pData(msms.dataset)$treat,mean)) )
apply(mSpC,2,summary)


###################################################
### code chunk number 4: PCAChunk
###################################################
snms <- substr(as.character(pData(msms.dataset)$treat),1,2)
snms <- paste(snms,as.integer(pData(msms.dataset)$batch),sep=".")
smpl.pca <- counts.pca(msms.dataset,snms=snms)$pca


###################################################
### code chunk number 5: msmsTests-Vignette2.Rnw:95-98
###################################################
par(mar=c(4,4,0.5,2)+0.1)
facs <- data.frame(batch=pData(msms.dataset)$batch)
counts.pca(msms.dataset,facs=facs,snms=snms)


###################################################
### code chunk number 6: SubsetChunk
###################################################
###  Subset and pre-process dataset
fl <- pData(msms.dataset)$batch=="2502"
e <- msms.dataset[,fl]
e <- pp.msms.data(e)
### Null and alternative model
null.f <- "y~1"
alt.f <- "y~treat"
### Normalizing condition
counts <- exprs(e)
div <- apply(counts,2,sum)
### Quasi-likelihood GLM
ql.res <- msms.glm.qlll(e,alt.f,null.f,div=div)
### Adjust p-values with FDR control.
adjp <- p.adjust(ql.res$p.value,method="BH")
###  Truth table
nh <- length(grep("HUMAN",featureNames(e)))
ny <- length(grep("HUMAN",featureNames(e),invert=TRUE))
tp <- length(grep("HUMAN",rownames(ql.res)[adjp<=0.05]))
fp <- sum(adjp<=0.05)-tp
(tt.ql1 <- data.frame(TP=tp,FP=fp,TN=ny-fp,FN=nh-tp))


###################################################
### code chunk number 7: SubsetFChunk
###################################################
###  Post-test filter
ql.tbl <- test.results(ql.res,e,pData(e)$treat,"U600","U200",div,
                    alpha=0.05,minSpC=2,minLFC=1,method="BH")$tres
ql.nms <- rownames(ql.tbl)[ql.tbl$DEP]
###  Truth table
ridx <- grep("HUMAN",ql.nms)
tp <- length(ridx)
fp <- length(ql.nms)-length(ridx)
(tt.ql11 <- data.frame(TP=tp,FP=fp,TN=ny-fp,FN=nh-tp))


###################################################
### code chunk number 8: GlobalChunk
###################################################
###  Pre-process dataset
gble <- pp.msms.data(msms.dataset)
### Null and alternative model
null.f <- "y~1"
alt.f <- "y~treat"
### Normalizing condition
div <- apply(exprs(gble),2,sum)
### Quasi-likelihood GLM
ql.res <- msms.glm.qlll(gble,alt.f,null.f,div=div)
### Adjust p-values with FDR control.
adjp <- p.adjust(ql.res$p.value,method="BH")
###  Truth table
nh <- length(grep("HUMAN",featureNames(gble)))
ny <- length(grep("HUMAN",featureNames(gble),invert=TRUE))
tp <- length(grep("HUMAN",rownames(ql.res)[adjp<=0.05]))
fp <- sum(adjp<=0.05)-tp
(tt.ql2 <- data.frame(TP=tp,FP=fp,TN=ny-fp,FN=nh-tp))


###################################################
### code chunk number 9: GlobalFChunk
###################################################
###  Post-test filter
ql.tbl <- test.results(ql.res,gble,pData(gble)$treat,"U600","U200",div,
                    alpha=0.05,minSpC=2,minLFC=1,method="BH")$tres
ql.nms <- rownames(ql.tbl)[ql.tbl$DEP]
###  Truth table
ridx <- grep("HUMAN",ql.nms)
tp <- length(ridx)
fp <- length(ql.nms)-length(ridx)
(tt.ql22 <- data.frame(TP=tp,FP=fp,TN=ny-fp,FN=nh-tp))


###################################################
### code chunk number 10: BlockChunk
###################################################
### Null and alternative model
null.f <- "y~batch"
alt.f <- "y~treat+batch"
### Quasi-likelihood GLM
ql.res <- msms.glm.qlll(gble,alt.f,null.f,div=div)
### Adjust p-values with FDR control.
adjp <- p.adjust(ql.res$p.value,method="BH")
###  Truth table
nh <- length(grep("HUMAN",featureNames(gble)))
ny <- length(grep("HUMAN",featureNames(gble),invert=TRUE))
tp <- length(grep("HUMAN",rownames(ql.res)[adjp<=0.05]))
fp <- sum(adjp<=0.05)-tp
(tt.ql3 <- data.frame(TP=tp,FP=fp,TN=ny-fp,FN=nh-tp))


###################################################
### code chunk number 11: BlockFChunk
###################################################
###  Post-test filter
ql.tbl <- test.results(ql.res,gble,pData(gble)$treat,"U600","U200",div,
                    alpha=0.05,minSpC=2,minLFC=1,method="BH")$tres
ql.nms <- rownames(ql.tbl)[ql.tbl$DEP]
###  Truth table
ridx <- grep("HUMAN",ql.nms)
tp <- length(ridx)
fp <- length(ql.nms)-length(ridx)
(tt.ql33 <- data.frame(TP=tp,FP=fp,TN=ny-fp,FN=nh-tp))


###################################################
### code chunk number 12: CompChunk
###################################################
gbl.tt <- rbind(tt.ql1,tt.ql11,tt.ql2,tt.ql22,tt.ql3,tt.ql33)
rownames(gbl.tt) <- c("subset","subset filtered",
                      "global","global filtered",
                      "blocking","blocking filtered")		  
library(xtable)
print(xtable(gbl.tt,align=c("r","r","r","r","r"),
             caption=c("Truth tables"),
             display=c("s","d","d","d","d")),
      type="latex",hline.after=c(-1,0,2,4,6),include.rownames=TRUE)


###################################################
### code chunk number 13: msmsTests-Vignette2.Rnw:256-263
###################################################
par(mar=c(5, 3, 2, 2) + 0.1)
par(cex.axis=0.8,cex.lab=0.8)
rownames(gbl.tt) <- c("subset","subset\nfiltered",
                      "global","global\nfiltered",
                      "blocking","blocking\nfiltered")		  
barplot(t(data.matrix(gbl.tt[,1:2])),beside=TRUE,las=2,space=c(0,0.25),
        legend.text=c("TP","FP"),args.legend=list(x="topleft",cex=0.8,bty="n"))