inst/doc/sSeq.R

In sSeq: Shrinkage estimation of dispersion in Negative Binomial models for RNA-seq experiments with small sample size

### R code from vignette source 'sSeq.Rnw'

###################################################
### code chunk number 1: hammer
###################################################
library(sSeq);
data(Hammer2months);
head(countsTable); 
conds.Hammer=c("A","A","B","B"); 


###################################################
### code chunk number 2: res1
###################################################
#exact test to get differential expressed genes
res1 = nbTestSH( countsTable, conds.Hammer, "A", "B");


###################################################
### code chunk number 3: seeResult
###################################################
head(res1);


###################################################
### code chunk number 4: AvsB_ASDplot
###################################################
disp1 <- nbTestSH( countsTable, conds.Hammer, "A", "B", SHonly=TRUE, plotASD=TRUE);


###################################################
### code chunk number 5: seeResult
###################################################
head(disp1);


###################################################
### code chunk number 6: dispPlot
###################################################
plotDispersion(disp1, legPos="bottomright")


###################################################
### code chunk number 7: variancePlot
###################################################
 rV = rowVars(countsTable);
 mu = rowMeans(countsTable); 
 SH.var =  (disp1$SH * mu^2 + mu)
 smoothScatter(log(rV)~log(mu), main="Variance Plot", ylab='log(variance)', 
   xlab='log(mean)', col=blues9[5], cex.axis=1.8)
 points( log(SH.var)~log(mu), col="black", pch=16)
 leg1 =  expression(paste("log(", hat("V")[g]^"sSeq", ")", sep=''));
 leg2 =  expression(paste("log(", hat("V")[g]^"MM", ")", sep=''));
 legend("bottomright", legend=c(leg1,leg2), col=c("black",blues9[5]), 
   pch=c(16, 1), cex=2)


###################################################
### code chunk number 8: ecdf
###################################################
#obtain the p-values for the comparison AvsA.
conds2.Hammer = c("A","B");
res1.2 =  nbTestSH( countsTable[,1:2], conds2.Hammer, "A", "B");
#draw the ECDF plot;
dd = data.frame(AvsA=res1.2$pval, AvsB=res1$pval);
ecdfAUC(dd, col.line=c("green", "red"), main = "ECDF, Hammer", drawRef = TRUE, rm1=TRUE)


###################################################
### code chunk number 9: MV_volcano
###################################################
drawMA_vol(countsTable, conds.Hammer, res1$pval, cutoff=0.05);


###################################################
### code chunk number 10: pairedDesign_data
###################################################
data(Tuch); 
head(countsTable); 
conds2 = c("normal", "normal", "normal", "tumor", "tumor", "tumor");
coLevels=data.frame(subjects=c(1,2,3,1,2,3));


###################################################
### code chunk number 11: pairedDesign
###################################################
res2 = nbTestSH(countsTable, conds2, "normal", "tumor", 
   coLevels= coLevels, pairedDesign=TRUE, pairedDesign.dispMethod="pooled");


###################################################
### code chunk number 12: seeResult
###################################################
head(res2)
countsTable[order(res2$pval),][1:25,]