inst/doc/msgbsR_Vignette.R

In msgbsR: msgbsR: methylation sensitive genotyping by sequencing (MS-GBS) R functions

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

###################################################
### code chunk number 1: load the example data
###################################################
library(msgbsR)
library(GenomicRanges)
library(SummarizedExperiment)

my_path <- system.file("extdata", package = "msgbsR")
se <- rawCounts(bamFilepath = my_path)
dim(assay(se))


###################################################
### code chunk number 2: insert metadata
###################################################
colData(se) <- DataFrame(Group = c(rep("Control", 3), rep("Experimental", 3)),
                         row.names = colnames(assay(se)))


###################################################
### code chunk number 3: extract cut sites
###################################################
cutSites <- rowRanges(se)

# # Adjust the cut sites to overlap recognition site on each strand
start(cutSites) <- ifelse(test = strand(cutSites) == '+',
                          yes = start(cutSites) - 1, no = start(cutSites) - 2)
end(cutSites) <- ifelse(test = strand(cutSites) == '+',
                          yes = end(cutSites) + 2, no = end(cutSites) + 1)


###################################################
### code chunk number 4: run checkCuts with BSgenome
###################################################
library(BSgenome.Rnorvegicus.UCSC.rn6)
correctCuts <- checkCuts(cutSites = cutSites, genome = "rn6", seq = "CCGG")


###################################################
### code chunk number 5: filter out incorrect cuts
###################################################
se <- subsetByOverlaps(se, correctCuts)
dim(assay(se))


###################################################
### code chunk number 6: plot counts per cut sites
###################################################
plotCounts(se = se, cateogory = "Group")


###################################################
### code chunk number 7: fig1
###################################################
plotCounts(se = se, cateogory = "Group")


###################################################
### code chunk number 8: differential methylation
###################################################
top <- diffMeth(se = se, cateogory = "Group",
                condition1 = "Control", condition2 = "Experimental",
                cpmThreshold = 1, thresholdSamples = 1)


###################################################
### code chunk number 9: chr20 length
###################################################
ratChr <- seqlengths(BSgenome.Rnorvegicus.UCSC.rn6)["chr20"]


###################################################
### code chunk number 10: top sites
###################################################
my_cuts <- GRanges(top$site[which(top$FDR < 0.05)])


###################################################
### code chunk number 11: circos plot
###################################################
plotCircos(cutSites = my_cuts, seqlengths = ratChr,
           cutSite.colour = "red", seqlengths.colour = "blue")


###################################################
### code chunk number 12: fig2
###################################################
plotCircos(cutSites = my_cuts, seqlengths = ratChr,
           cutSite.colour = "red", seqlengths.colour = "blue")


###################################################
### code chunk number 13: annotation
###################################################
sessionInfo()