### R code from vignette source 'DSS.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex(use.unsrturl=FALSE) ################################################### ### code chunk number 2: DSS.Rnw:119-122 (eval = FALSE) ################################################### ## library(GenomicFeatures) ## txdb = makeTranscriptDbFromUCSC(genom="hg19",tablename="refGene") ## genes = genes(txdb) ################################################### ### code chunk number 3: DSS.Rnw:126-127 (eval = FALSE) ################################################### ## bam=scanBam("data.bam") ################################################### ### code chunk number 4: DSS.Rnw:130-131 (eval = FALSE) ################################################### ## IRange.reads=GRanges(seqnames=Rle(bam$rname), ranges=IRanges(bam$pos, width=bam$qwidth)) ################################################### ### code chunk number 5: DSS.Rnw:135-136 (eval = FALSE) ################################################### ## counts = countOverlaps(genes, IRange.reads) ################################################### ### code chunk number 6: DSS.Rnw:170-179 ################################################### 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 7: DSS.Rnw:182-183 ################################################### seqData=estNormFactors(seqData) ################################################### ### code chunk number 8: DSS.Rnw:186-187 ################################################### seqData=estDispersion(seqData) ################################################### ### code chunk number 9: DSS.Rnw:191-193 ################################################### result=waldTest(seqData, 0, 1) head(result,5) ################################################### ### code chunk number 10: DSS.Rnw:201-205 ################################################### counts = matrix(rpois(600, 10), ncol=6) designs = c(0,0,0,1,1,1) result = DSS.DE(counts, designs) head(result) ################################################### ### code chunk number 11: DSS.Rnw:220-225 ################################################### 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 12: DSS.Rnw:229-232 ################################################### seqData=newSeqCountSet(counts, as.data.frame(X)) seqData=estNormFactors(seqData) seqData=estDispersion(seqData) ################################################### ### code chunk number 13: DSS.Rnw:238-239 ################################################### fit.edgeR <- glmFit(counts, X, dispersion=dispersion(seqData)) ################################################### ### code chunk number 14: DSS.Rnw:243-245 ################################################### lrt.edgeR <- glmLRT(glmfit=fit.edgeR, coef=2) head(lrt.edgeR$table) ################################################### ### code chunk number 15: DSS.Rnw:338-348 ################################################### 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) BSobj <- makeBSseqData( list(dat1.1, dat1.2, dat2.1, dat2.2), c("C1","C2", "N1", "N2") )[1:1000,] BSobj ################################################### ### code chunk number 16: DSS.Rnw:361-363 ################################################### dmlTest <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2")) head(dmlTest) ################################################### ### code chunk number 17: DSS.Rnw:367-368 ################################################### dmlTest.sm <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"), smoothing=TRUE) ################################################### ### code chunk number 18: DSS.Rnw:386-388 ################################################### dmls <- callDML(dmlTest, p.threshold=0.001) head(dmls) ################################################### ### code chunk number 19: DSS.Rnw:394-396 ################################################### dmls2 <- callDML(dmlTest, delta=0.1, p.threshold=0.001) head(dmls2) ################################################### ### code chunk number 20: DSS.Rnw:410-412 ################################################### dmrs <- callDMR(dmlTest, p.threshold=0.01) head(dmrs) ################################################### ### code chunk number 21: DSS.Rnw:419-421 ################################################### dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05) head(dmrs2) ################################################### ### code chunk number 22: DSS.Rnw:437-438 (eval = FALSE) ################################################### ## showOneDMR(dmrs[1,], BSobj) ################################################### ### code chunk number 23: DSS.Rnw:480-484 ################################################### Strain = rep(c("A", "B", "C"), 4) Sex = rep(c("M", "F"), each=6) design = data.frame(Strain,Sex) design ################################################### ### code chunk number 24: DSS.Rnw:489-491 ################################################### X = model.matrix(~Strain+ Sex, design) X ################################################### ### code chunk number 25: DSS.Rnw:509-511 ################################################### L = cbind(c(0,1,0,0),c(0,0,1,0)) L ################################################### ### code chunk number 26: DSS.Rnw:517-518 ################################################### matrix(c(0,1,-1,0), ncol=1) ################################################### ### code chunk number 27: DSS.Rnw:532-535 ################################################### data(RRBS) RRBS design ################################################### ### code chunk number 28: DSS.Rnw:540-541 ################################################### DMLfit = DMLfit.multiFactor(RRBS, design=design, formula=~case+cell+case:cell) ################################################### ### code chunk number 29: DSS.Rnw:551-552 ################################################### DMLtest.cell = DMLtest.multiFactor(DMLfit, coef=3) ################################################### ### code chunk number 30: DSS.Rnw:559-560 ################################################### colnames(DMLfit$X) ################################################### ### code chunk number 31: DSS.Rnw:566-567 ################################################### DMLtest.cell = DMLtest.multiFactor(DMLfit, coef="cellrN") ################################################### ### code chunk number 32: DSS.Rnw:575-577 ################################################### ix=sort(DMLtest.cell[,"pvals"], index.return=TRUE)$ix head(DMLtest.cell[ix,]) ################################################### ### code chunk number 33: DSS.Rnw:582-585 ################################################### par(mfrow=c(1,2)) hist(DMLtest.cell$stat, 50, main="test statistics", xlab="") hist(DMLtest.cell$pvals, 50, main="P values", xlab="") ################################################### ### code chunk number 34: DSS.Rnw:589-590 ################################################### callDMR(DMLtest.cell, p.threshold=0.05) ################################################### ### code chunk number 35: DSS.Rnw:604-613 ################################################### ## fit a model with additive effect only DMLfit = DMLfit.multiFactor(RRBS, design, ~case+cell) ## test case effect test1 = DMLtest.multiFactor(DMLfit, coef=2) test2 = DMLtest.multiFactor(DMLfit, coef="caseSLE") test3 = DMLtest.multiFactor(DMLfit, term="case") Contrast = matrix(c(0,1,0), ncol=1) test4 = DMLtest.multiFactor(DMLfit, Contrast=Contrast) cor(cbind(test1$pval, test2$pval, test3$pval, test4$pval)) ################################################### ### code chunk number 36: DSS.Rnw:625-626 ################################################### sessionInfo()