### R code from vignette source 'vignettes/DSS/inst/doc/DSS.Rnw' ################################################### ### code chunk number 1: DSS.Rnw:132-141 ################################################### library(DSS) counts1=matrix(rnbinom(300, mu=10, size=10), ncol=3) counts2=matrix(rnbinom(300, mu=50, size=10), ncol=3) X1=cbind(counts1, counts2) ## these are 100 DE genes X2=matrix(rnbinom(11400, mu=10, size=10), ncol=6) X=rbind(X1,X2) designs=c(0,0,0,1,1,1) seqData=newSeqCountSet(X, designs) seqData ################################################### ### code chunk number 2: DSS.Rnw:144-145 ################################################### seqData=estNormFactors(seqData) ################################################### ### code chunk number 3: DSS.Rnw:148-149 ################################################### seqData=estDispersion(seqData) ################################################### ### code chunk number 4: DSS.Rnw:153-155 ################################################### result=waldTest(seqData, 0, 1) head(result,5) ################################################### ### code chunk number 5: DSS.Rnw:170-175 ################################################### library(DSS) library(edgeR) counts=matrix(rpois(800, 10), ncol=8) design=data.frame(gender=c(rep("M",4), rep("F",4)), strain=rep(c("WT", "Mutant"),4)) X=model.matrix(~gender+strain, data=design) ################################################### ### code chunk number 6: DSS.Rnw:179-182 ################################################### seqData=newSeqCountSet(counts, as.data.frame(X)) seqData=estNormFactors(seqData) seqData=estDispersion(seqData) ################################################### ### code chunk number 7: DSS.Rnw:187-189 ################################################### fit.edgeR <- glmFit(counts, X, lib.size=normalizationFactor(seqData), dispersion=dispersion(seqData)) ################################################### ### code chunk number 8: DSS.Rnw:193-195 ################################################### lrt.edgeR <- glmLRT(fit.edgeR, coef=2) head(lrt.edgeR$table) ################################################### ### code chunk number 9: DSS.Rnw:220-230 ################################################### library(DSS) require(bsseq) path <- file.path(system.file(package="DSS"), "extdata") dat1.1 <- read.table(file.path(path, "cond1_1.txt"), header=TRUE) dat1.2 <- read.table(file.path(path, "cond1_2.txt"), header=TRUE) dat2.1 <- read.table(file.path(path, "cond2_1.txt"), header=TRUE) dat2.2 <- read.table(file.path(path, "cond2_2.txt"), header=TRUE) BS1 <- makeBSseqData( list(dat1.1, dat1.2), paste("cond1",1:2,sep=".") ) BS2 <- makeBSseqData( list(dat2.1, dat2.2), paste("cond2",1:2,sep=".") ) BS1 ################################################### ### code chunk number 10: DSS.Rnw:234-236 ################################################### dmls <- callDML(BS1, BS2) head(dmls) ################################################### ### code chunk number 11: DSS.Rnw:240-242 ################################################### dmrs <- callDMR(dmls, p.threshold=0.001) head(dmrs) ################################################### ### code chunk number 12: DSS.Rnw:262-263 ################################################### sessionInfo()