inst/doc/RJMCMCNucleosomes.R

In RJMCMCNucleosomes: Bayesian hierarchical model for genome-wide nucleosome positioning with high-throughput short-read data (MNase-Seq)

## ----loadingPackage, warning=FALSE, message=FALSE-----------------------------
library(RJMCMCNucleosomes)

## ----createSample, collapse=TRUE, message=FALSE-------------------------------
## Load nucleoSim package
library(nucleoSim)

val.num       <- 50     ### Number of well-positioned nucleosomes
val.del       <- 10     ### Number of well-positioned nucleosomes to delete
val.var       <- 30     ### variance associated to well-positioned nucleosomes
val.fuz       <- 10     ### Number of fuzzy nucleosomes
val.fuz.var   <- 50     ### variance associated to fuzzy nucleosomes
val.max.cover <- 70     ### Maximum coverage for one nucleosome
val.nuc.len   <- 147    ### Distance between nucleosomes
val.len.var   <- 10     ### Variance associated to the length of the reads
val.lin.len   <- 20     ### The length of the DNA linker
val.rnd.seed  <- 100    ### Set seed when result needs to be reproducible
val.offset    <- 10000  ### The number of bases used to offset 
                        ### all nucleosomes and reads

## Create sample using a Normal distribution
sample <- nucleoSim::syntheticNucReadsFromDist(wp.num=val.num,
                                    wp.del=val.del, wp.var=val.var,
                                    fuz.num=val.del, fuz.var=val.fuz.var,
                                    max.cover=val.max.cover, 
                                    nuc.len=val.nuc.len,
                                    len.var=val.len.var, 
                                    lin.len=val.lin.len,
                                    rnd.seed=val.rnd.seed,
                                    distr="Normal", offset=val.offset)

## Create visual representation of the synthetic nucleosome sample
plot(sample)

## ----segment01, warning=FALSE, collapse=TRUE, message=FALSE-------------------
## Load needed packages
library(GenomicRanges)

## Transform sample dataset into GRanges object
sampleGRanges <- GRanges(seqnames = sample$dataIP$chr, 
                        ranges = IRanges(start = sample$dataIP$start, 
                                        end = sample$dataIP$end), 
                        strand = sample$dataIP$strand)

## Segment sample into candidate regions
sampleSegmented <- segmentation(reads = sampleGRanges, zeta = 147, 
                                delta = 40, maxLength = 1000)

## Number of segments created
length(sampleSegmented)

## ----runRJMCMC01 , warning=FALSE, collapse=TRUE-------------------------------
## Extract the first segment 
segment01 <- sampleSegmented[[1]]

## Run RJMCMC analysis
## A higher number of iterations is recommanded for real analysis
resultSegment01 <- rjmcmc(reads  = segment01, nbrIterations = 100000, 
                            lambda = 3, kMax = 30,
                            minInterval = 100, maxInterval = 200, 
                            minReads = 5, vSeed = 1921)

## Print the predicted nucleosomes for the first segment
resultSegment01

## ----regroup01, warning=FALSE, collapse=TRUE, message=FALSE-------------------
## The directory containing the results of all segments
## On RDS file has been created for each segment
directory <- system.file("extdata", "demo_vignette", 
                            package = "RJMCMCNucleosomes")

## Merging predicted nucleosomes from all segments
resultsAllMerged <- mergeAllRDSFilesFromDirectory(directory)

resultsAllMerged

## ----postProcess01, collapse=TRUE, message=FALSE------------------------------
## Split reads from the initial sample data into forward and reverse subsets
reads <- GRanges(sample$dataIP)

## Number of nucleosomes before the post-treatment
resultsAllMerged$k

## Use the post-treatment function
resultsPostTreatment <- postTreatment(reads = reads,
                            resultRJMCMC = resultsAllMerged,
                            extendingSize = 15,
                            chrLength = max(start(reads), end(reads)) + 1000)


## Number of nucleosomes after the post-treatment
length(resultsPostTreatment)

## ----visualisation, collapse=TRUE, message=FALSE, fig.height=6.5, fig.width=8----
## Extract reads to create a GRanges
reads <-GRanges(sample$dataIP)

## Merge predictions from before post-treatment and after post-treatment in 
## a list so that both results will be shown in graph
# resultsBeforeAndAfter <- list(Sample = c(sample$wp$nucleopos,
#                                             sample$fuz$nucleopos),
#                                BeforePostTreatment = resultsAllMerged$mu,
#                                AfterPostTreatment = resultsPostTreatment)
resultsBeforeAndAfter <- GRangesList(Sample = GRanges(
        rep("chr_SYNTHETIC", 
            length(c(sample$wp$nucleopos,sample$fuz$nucleopos))), 
        ranges = IRanges(start=c(sample$wp$nucleopos,sample$fuz$nucleopos),
                        end=c(sample$wp$nucleopos,sample$fuz$nucleopos)), 
        strand=rep("*", length(c(sample$wp$nucleopos,
                                sample$fuz$nucleopos)))),
                        BeforePostTreatment = resultsAllMerged$mu,
                        AfterPostTreatment = resultsPostTreatment)
## Create graph using nucleosome positions and reads
## The plot will shows : 
## 1. nucleosomes from the sample
## 2. nucleosomes detected by rjmcmc() function
## 3. nucleosomes obtained after post-treament 
plotNucleosomes(nucleosomePositions = resultsBeforeAndAfter, reads = reads, 
                    names=c("Sample", "RJMCMC", "After Post-Treatment"))

## ----rjmcmcCHR, collapse=TRUE, message=FALSE, eval=FALSE----------------------
#  ## Load synthetic dataset of reads
#  data(syntheticNucleosomeReads)
#  
#  ## Number of reads in the dataset
#  nrow(syntheticNucleosomeReads$dataIP)
#  
#  ## Use the dataset to create a GRanges object
#  sampleGRanges <- GRanges(syntheticNucleosomeReads$dataIP)
#  
#  ## All reads are related to one chromosome called "chr_SYNTHETIC"
#  seqnames(sampleGRanges)
#  
#  ## Run RJMCMC on all reads
#  result <- rjmcmcCHR(reads = sampleGRanges,
#                  seqName = "chr_SYNTHETIC", dirOut = "testRJMCMCCHR",
#                  zeta = 147, delta=50, maxLength=1200,
#                  nbrIterations = 500, lambda = 3, kMax = 30,
#                  minInterval = 146, maxInterval = 292, minReads = 5,
#                  vSeed = 10113, nbCores = 2, saveAsRDS = FALSE,
#                  saveSEG = FALSE)
#  
#  result

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()