BiSeq.R
In BiSeq: Processing and analyzing bisulfite sequencing data

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

###################################################
### code chunk number 1: BiSeq.Rnw:35-37
###################################################
options(width=60)
options(continue=" ")


###################################################
### code chunk number 2: preliminaries
###################################################
library(BiSeq)


###################################################
### code chunk number 3: read (eval = FALSE)
###################################################
## readBismark(files, colData)


###################################################
### code chunk number 4: BiSeq.Rnw:81-93
###################################################
metadata <- list(Sequencer = "Instrument", Year = "2013")
rowRanges <- GRanges(seqnames = "chr1",
                  ranges = IRanges(start = c(1,2,3), end = c(1,2,3)))
colData <- DataFrame(group = c("cancer", "control"),
                     row.names = c("sample_1", "sample_2"))
totalReads <- matrix(c(rep(10L, 3), rep(5L, 3)), ncol = 2)
methReads <- matrix(c(rep(5L, 3), rep(5L, 3)), ncol = 2)
BSraw(metadata = metadata,
      rowRanges = rowRanges,
      colData = colData,
      totalReads = totalReads,
      methReads = methReads)


###################################################
### code chunk number 5: BiSeq.Rnw:98-100
###################################################
data(rrbs)
rrbs


###################################################
### code chunk number 6: BiSeq.Rnw:103-104
###################################################
colData(rrbs)


###################################################
### code chunk number 7: BiSeq.Rnw:107-108
###################################################
head(rowRanges(rrbs))


###################################################
### code chunk number 8: BiSeq.Rnw:111-112
###################################################
head(totalReads(rrbs))


###################################################
### code chunk number 9: BiSeq.Rnw:115-116
###################################################
head(methReads(rrbs))


###################################################
### code chunk number 10: BiSeq.Rnw:132-137
###################################################
methLevel <- matrix(c(rep(0.5, 3), rep(1, 3)), ncol = 2)
BSrel(metadata = metadata,
      rowRanges = rowRanges,
      colData = colData,
      methLevel = methLevel)


###################################################
### code chunk number 11: BiSeq.Rnw:141-143
###################################################
rrbs.rel <- rawToRel(rrbs)
rrbs.rel


###################################################
### code chunk number 12: BiSeq.Rnw:146-147
###################################################
head(methLevel(rrbs.rel))


###################################################
### code chunk number 13: BiSeq.Rnw:152-154
###################################################
dim(rrbs)
colnames(rrbs)


###################################################
### code chunk number 14: BiSeq.Rnw:157-160
###################################################
rrbs[,"APL2"]
ind.chr1 <- which(seqnames(rrbs) == "chr1")
rrbs[ind.chr1,]


###################################################
### code chunk number 15: BiSeq.Rnw:163-166
###################################################
region <- GRanges(seqnames="chr1",
                  ranges=IRanges(start = 875200,
                                 end = 875500))


###################################################
### code chunk number 16: BiSeq.Rnw:168-170
###################################################
findOverlaps(rrbs, region)
subsetByOverlaps(rrbs, region)


###################################################
### code chunk number 17: BiSeq.Rnw:173-174
###################################################
sort(rrbs)


###################################################
### code chunk number 18: BiSeq.Rnw:177-180
###################################################
combine(rrbs[1:10,1:2], rrbs[1:1000, 3:10])
split(rowRanges(rrbs),
      f = as.factor(as.character(seqnames(rrbs))))


###################################################
### code chunk number 19: BiSeq.Rnw:186-187
###################################################
covStatistics(rrbs)


###################################################
### code chunk number 20: BiSeq.Rnw:191-192
###################################################
covBoxplots(rrbs, col = "cornflowerblue", las = 2)


###################################################
### code chunk number 21: BiSeq.Rnw:218-224
###################################################
rrbs.small <- rrbs[1:1000,]
rrbs.clust.unlim <- clusterSites(object = rrbs.small,
                                 groups = colData(rrbs)$group,
                                 perc.samples = 4/5,
                                 min.sites = 20,
                                 max.dist = 100)


###################################################
### code chunk number 22: BiSeq.Rnw:227-228
###################################################
head(rowRanges(rrbs.clust.unlim))


###################################################
### code chunk number 23: BiSeq.Rnw:231-232
###################################################
clusterSitesToGR(rrbs.clust.unlim)


###################################################
### code chunk number 24: BiSeq.Rnw:237-241
###################################################
ind.cov <- totalReads(rrbs.clust.unlim) > 0
quant <- quantile(totalReads(rrbs.clust.unlim)[ind.cov], 0.9)
quant
rrbs.clust.lim <- limitCov(rrbs.clust.unlim, maxCov = quant)


###################################################
### code chunk number 25: BiSeq.Rnw:245-246
###################################################
covBoxplots(rrbs.clust.lim, col = "cornflowerblue", las = 2)


###################################################
### code chunk number 26: BiSeq.Rnw:252-253
###################################################
predictedMeth <- predictMeth(object = rrbs.clust.lim)


###################################################
### code chunk number 27: BiSeq.Rnw:256-257
###################################################
predictedMeth


###################################################
### code chunk number 28: BiSeq.Rnw:262-268
###################################################
plotMeth(object.raw = rrbs[,6],
         object.rel = predictedMeth[,6],
         region = region,
         lwd.lines = 2,
         col.points = "blue",
         cex = 1.5)


###################################################
### code chunk number 29: BiSeq.Rnw:279-288
###################################################
cancer <- predictedMeth[, colData(predictedMeth)$group == "APL"]
control <- predictedMeth[, colData(predictedMeth)$group == "control"]
mean.cancer <- rowMeans(methLevel(cancer))
mean.control <- rowMeans(methLevel(control))
plot(mean.control,
     mean.cancer,
     col = "blue",
     xlab = "Methylation in controls",
     ylab = "Methylation in APLs")


###################################################
### code chunk number 30: BiSeq.Rnw:294-299
###################################################
## To shorten the run time set mc.cores, if possible!
betaResults <- betaRegression(formula = ~group,
                              link = "probit",
                              object = predictedMeth,
                              type = "BR")


###################################################
### code chunk number 31: BiSeq.Rnw:301-303
###################################################
## OR:
data(betaResults)


###################################################
### code chunk number 32: BiSeq.Rnw:306-307
###################################################
head(betaResults)


###################################################
### code chunk number 33: BiSeq.Rnw:314-323
###################################################
## Both resampled groups should have the same number of
## cancer and control samples:
predictedMethNull <- predictedMeth[,c(1:4, 6:9)]
colData(predictedMethNull)$group.null <- rep(c(1,2), 4)
## To shorten the run time, please set mc.cores, if possible!
betaResultsNull <- betaRegression(formula = ~group.null,
                                  link = "probit",
                                  object = predictedMethNull,
                                  type="BR")


###################################################
### code chunk number 34: BiSeq.Rnw:325-327
###################################################
## OR:
data(betaResultsNull)


###################################################
### code chunk number 35: BiSeq.Rnw:330-331
###################################################
vario <- makeVariogram(betaResultsNull)


###################################################
### code chunk number 36: BiSeq.Rnw:333-335
###################################################
## OR:
data(vario)


###################################################
### code chunk number 37: BiSeq.Rnw:340-345
###################################################
plot(vario$variogram$v)
vario.sm <- smoothVariogram(vario, sill = 0.9)
lines(vario.sm$variogram[,c("h", "v.sm")],
      col = "red", lwd = 1.5)
grid()


###################################################
### code chunk number 38: BiSeq.Rnw:351-354
###################################################
names(vario.sm)
head(vario.sm$variogram)
head(vario.sm$pValsList[[1]])


###################################################
### code chunk number 39: BiSeq.Rnw:357-362
###################################################
## auxiliary object to get the pValsList for the test
## results of interest:
vario.aux <- makeVariogram(betaResults, make.variogram=FALSE)
vario.sm$pValsList <- vario.aux$pValsList
head(vario.sm$pValsList[[1]])


###################################################
### code chunk number 40: BiSeq.Rnw:365-366
###################################################
locCor <- estLocCor(vario.sm)


###################################################
### code chunk number 41: BiSeq.Rnw:369-372
###################################################
clusters.rej <- testClusters(locCor,
                             FDR.cluster = 0.1)
clusters.rej$clusters.reject


###################################################
### code chunk number 42: BiSeq.Rnw:379-382
###################################################
clusters.trimmed <- trimClusters(clusters.rej,
                                 FDR.loc = 0.05)
head(clusters.trimmed)


###################################################
### code chunk number 43: BiSeq.Rnw:389-393
###################################################
DMRs <- findDMRs(clusters.trimmed,
                 max.dist = 100,
                 diff.dir = TRUE)
DMRs


###################################################
### code chunk number 44: BiSeq.Rnw:399-404
###################################################
DMRs.2 <- compareTwoSamples(object = predictedMeth,
                            sample1 = "APL1",
                            sample2 = "APL10961",
                            minDiff = 0.3,
                            max.dist = 100)


###################################################
### code chunk number 45: BiSeq.Rnw:407-408
###################################################
sum(overlapsAny(DMRs.2,DMRs))


###################################################
### code chunk number 46: <
###################################################
data(promoters)
data(rrbs)
rrbs.red <- subsetByOverlaps(rrbs, promoters)
ov <- findOverlaps(rrbs.red, promoters)
rowRanges(rrbs.red)$cluster.id[queryHits(ov)] <- promoters$acc_no[subjectHits(ov)]
head(rowRanges(rrbs.red))


###################################################
### code chunk number 47: <
###################################################
data(promoters)
data(rrbs)
rrbs <- rawToRel(rrbs)
promoters <- promoters[overlapsAny(promoters, rrbs)]
gt <- globalTest(group~1,
                 rrbs,
                 subsets = promoters)
head(gt)


###################################################
### code chunk number 48: BiSeq.Rnw:473-481
###################################################
rowCols <- c("magenta", "blue")[as.numeric(colData(predictedMeth)$group)]
plotMethMap(predictedMeth,
            region = DMRs[3],
            groups = colData(predictedMeth)$group,
            intervals = FALSE,
            zlim = c(0,1),
            RowSideColors = rowCols,
            labCol = "", margins = c(0, 6))


###################################################
### code chunk number 49: BiSeq.Rnw:489-499
###################################################
plotSmoothMeth(object.rel = predictedMeth,
               region = DMRs[3],
               groups = colData(predictedMeth)$group,
               group.average = FALSE,
               col = c("magenta", "blue"),
               lwd = 1.5)
legend("topright",
       legend=levels(colData(predictedMeth)$group),
       col=c("magenta", "blue"),
       lty=1, lwd = 1.5)


###################################################
### code chunk number 50: BiSeq.Rnw:506-513
###################################################
data(promoters)
head(promoters)
DMRs.anno <- annotateGRanges(object = DMRs,
                             regions = promoters,
                             name = 'Promoter',
                             regionInfo = 'acc_no')
DMRs.anno


###################################################
### code chunk number 51: BiSeq.Rnw:519-529
###################################################
plotBindingSites(object = rrbs,
                 regions = promoters,
                 width = 4000,
                 group = colData(rrbs)$group,
                 col = c("magenta", "blue"),
                 lwd = 1.5)
legend("top",
       legend=levels(colData(rrbs)$group),
       col=c("magenta", "blue"),
       lty=1, lwd = 1.5)


###################################################
### code chunk number 52: BiSeq.Rnw:535-545 (eval = FALSE)
###################################################
## track.names <- paste(colData(rrbs)$group,
##                      "_",
##                      gsub("APL", "", colnames(rrbs)),
##                      sep="")
## writeBED(object = rrbs,
##          name = track.names,
##          file = paste(colnames(rrbs), ".bed", sep = ""))
## writeBED(object = predictedMeth,
##          name = track.names,
##          file = paste(colnames(predictedMeth), ".bed", sep = ""))