inst/doc/BiologicalApplication.R

In MIRA: Methylation-Based Inference of Regulatory Activity

## ---- eval=TRUE, message=FALSE, warning=FALSE---------------------------------
library(MIRA)
library(data.table) # for the functions: fread, setkey, merge
library(GenomicRanges) # for the functions: GRanges, resize
library(ggplot2) # for the function: ylab
exampleAnnoDT2 <- fread(system.file("extdata", "exampleAnnoDT2.txt", 
                                   package="MIRA")) 

## ---- eval=FALSE--------------------------------------------------------------
#  # 12 Ewing samples: T1-T12
#  pathToData <- "/path/to/MIRA_Ewing_Vignette_Files/"
#  ewingFiles <- paste0(pathToData, "EWS_T", 1:12, ".bed")
#  # 12 muscle related samples, 3 of each type
#  muscleFiles <- c("Hsmm_", "Hsmmt_", "Hsmmfshd_","Hsmmtubefshd_")
#  muscleFiles <- paste0(pathToData, muscleFiles, rep(1:3, each=4), ".bed")
#  RRBSFiles <- c(ewingFiles, muscleFiles)
#  regionSetFiles <- paste0(pathToData, c("sknmc_specific.bed", "muscle_specific.bed"))

## ---- eval=FALSE--------------------------------------------------------------
#  BSDTList <- lapply(X=RRBSFiles, FUN=BSreadBiSeq)
#  regionSets <- lapply(X=regionSetFiles, FUN=fread)

## ---- eval=FALSE--------------------------------------------------------------
#  names(BSDTList) <- tools::file_path_sans_ext(basename(RRBSFiles))
#  regionSets <- lapply(X=regionSets,
#                      FUN=function(x) setNames(x, c("chr", "start", "end")))
#  regionSets <- lapply(X=regionSets, FUN=
#                          function(x) GRanges(seqnames=x$chr,
#                                                ranges=IRanges(x$start, x$end)))

## ---- eval=FALSE--------------------------------------------------------------
#  regionSets[[1]] <- resize(regionSets[[1]], 5000, fix="center")
#  regionSets[[2]] <- resize(regionSets[[2]], 5000, fix="center")
#  names(regionSets) <- c("Ewing_Specific", "Muscle_Specific")

## ---- eval=FALSE--------------------------------------------------------------
#  subBSDTList <- BSDTList[c(1, 5, 9, 13, 17, 21)]
#  bigBin <- lapply(X=subBSDTList, FUN=aggregateMethyl, GRList=regionSets,
#                  binNum=21, minBaseCovPerBin=0)

## ---- eval=FALSE--------------------------------------------------------------
#  bigBinDT1 <- rbindNamedList(bigBin, newColName = "sampleName")

## ---- eval=TRUE---------------------------------------------------------------
# from data included with MIRA package
data(bigBinDT1)

## ---- eval=TRUE---------------------------------------------------------------
setkey(exampleAnnoDT2, sampleName)
setkey(bigBinDT1, sampleName)
bigBinDT1 <- merge(bigBinDT1, exampleAnnoDT2, all.x=TRUE)

## ---- eval=TRUE---------------------------------------------------------------
plotMIRAProfiles(binnedRegDT=bigBinDT1)

## ---- eval=FALSE--------------------------------------------------------------
#  regionSets[[1]] <- resize(regionSets[[1]], 4000, fix="center") # Ewing
#  regionSets[[2]] <- resize(regionSets[[2]], 1750, fix="center") # muscle

## ---- eval=FALSE--------------------------------------------------------------
#  bigBin <- lapply(X=BSDTList, FUN=aggregateMethyl, GRList=regionSets,
#                  binNum=21, minBaseCovPerBin=0)
#  bigBinDT2 <- rbindNamedList(bigBin)

## ---- eval=TRUE---------------------------------------------------------------
# from data included with MIRA package
data(bigBinDT2)

## ---- eval=TRUE---------------------------------------------------------------
# data.table functions are used to add info from annotation object to bigBinDT2
setkey(exampleAnnoDT2, sampleName)
setkey(bigBinDT2, sampleName)
bigBinDT2 <- merge(bigBinDT2, exampleAnnoDT2, all.x=TRUE)

## ---- eval=TRUE---------------------------------------------------------------
plotMIRAProfiles(binnedRegDT=bigBinDT2)

## ---- eval=TRUE, results="hide"-----------------------------------------------
bigBinDT2[, methylProp := methylProp - min(methylProp) + .05, by=.(featureID, sampleName)]

## ---- eval=TRUE---------------------------------------------------------------
normPlot <- plotMIRAProfiles(binnedRegDT=bigBinDT2)
normPlot + ylab("Normalized DNA Methylation (%)")

## ---- eval=TRUE---------------------------------------------------------------
sampleScores <- calcMIRAScore(bigBinDT2,
                        shoulderShift="auto",
                        regionSetIDColName="featureID",
                        sampleIDColName="sampleName")
head(sampleScores)

## ---- eval=TRUE---------------------------------------------------------------
setkey(exampleAnnoDT2, sampleName)
setkey(sampleScores, sampleName)
sampleScores <- merge(sampleScores, exampleAnnoDT2, all.x=TRUE)
plotMIRAScores(sampleScores)

## -----------------------------------------------------------------------------
EwingSampleEwingRegions <- sampleScores[sampleScores$sampleType == "Ewing" &
                                           sampleScores$featureID ==
                                           "Ewing_Specific", ]
myoSampleEwingRegions <- sampleScores[sampleScores$sampleType == "Muscle-related" &
                                           sampleScores$featureID ==
                                           "Ewing_Specific", ]
wilcox.test(EwingSampleEwingRegions$score, myoSampleEwingRegions$score)

## -----------------------------------------------------------------------------
EwingSampleMyoRegions <- sampleScores[sampleScores$sampleType == "Ewing" &
                                           sampleScores$featureID ==
                                           "Muscle_Specific", ]
myoSampleMyoRegions <- sampleScores[sampleScores$sampleType == "Muscle-related" &
                                           sampleScores$featureID ==
                                           "Muscle_Specific", ]
wilcox.test(EwingSampleMyoRegions$score, myoSampleMyoRegions$score)