### 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-240
###################################################
fit.edgeR <- glmFit(counts, X, lib.size=normalizationFactor(seqData),
             dispersion=dispersion(seqData))


###################################################
### code chunk number 14: DSS.Rnw:244-246
###################################################
lrt.edgeR <- glmLRT(glmfit=fit.edgeR, coef=2)
head(lrt.edgeR$table)


###################################################
### code chunk number 15: DSS.Rnw:339-349
###################################################
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:10000,]
BSobj


###################################################
### code chunk number 16: DSS.Rnw:362-364
###################################################
dmlTest <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"))
head(dmlTest)


###################################################
### code chunk number 17: DSS.Rnw:368-369
###################################################
dmlTest.sm <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"), smoothing=TRUE)


###################################################
### code chunk number 18: DSS.Rnw:387-389
###################################################
  dmls <- callDML(dmlTest, p.threshold=0.001)
  head(dmls)


###################################################
### code chunk number 19: DSS.Rnw:395-397
###################################################
  dmls2 <- callDML(dmlTest, delta=0.1, p.threshold=0.001)
  head(dmls2)


###################################################
### code chunk number 20: DSS.Rnw:411-413
###################################################
dmrs <- callDMR(dmlTest, p.threshold=0.01)
head(dmrs)


###################################################
### code chunk number 21: DSS.Rnw:420-422
###################################################
  dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05)
  head(dmrs2)


###################################################
### code chunk number 22: DSS.Rnw:438-439 (eval = FALSE)
###################################################
##   showOneDMR(dmrs[1,], BSobj)


###################################################
### code chunk number 23: DSS.Rnw:469-472
###################################################
data(RRBS)
RRBS
design


###################################################
### code chunk number 24: DSS.Rnw:477-478
###################################################
DMLfit = DMLfit.multiFactor(RRBS, design=design, formula=~case+cell+case:cell)


###################################################
### code chunk number 25: DSS.Rnw:488-489
###################################################
DMLtest.cell = DMLtest.multiFactor(DMLfit, coef=3)


###################################################
### code chunk number 26: DSS.Rnw:497-499
###################################################
ix=sort(DMLtest.cell[,"pvals"], index.return=TRUE)$ix
head(DMLtest.cell[ix,])


###################################################
### code chunk number 27: DSS.Rnw:504-507
###################################################
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 28: DSS.Rnw:519-520
###################################################
sessionInfo()