## ----style, echo = FALSE, results = 'asis'------------------------------- BiocStyle::markdown() ## ---- message=FALSE, warning=FALSE--------------------------------------- library(minfi) library(minfiData) library(sva) library(Gviz) library(FlowSorted.Blood.450k) library(shinyMethyl) library(shinyMethylData) library(tutorial.450k) ## ----eval=FALSE---------------------------------------------------------- # # WARNING: No need to run this chunk if you're in the workshop at BioC 2016 and # # using the Amazon Machine Image # source("http://www.bioconductor.org/biocLite.R") # useDevel(TRUE) # # NOTE: You may not need `type = "source"` # biocLite(c("minfi", "minfiData", "sva", "Gviz", "FlowSorted.Blood.450k", # "shinyMethyl", "shinyMethylData"), # type = "source") ## ----eval=FALSE---------------------------------------------------------- # # NOTE: You'll also need the devtools package to install the tutorial from # # GitHub # biocLite("devtools") # biocLite("hansenlab/tutorial.450k") ## ------------------------------------------------------------------------ baseDir <- system.file("extdata", package="minfiData") targets <- read.metharray.sheet(baseDir) targets RGSet <- read.metharray.exp(targets = targets) ## ------------------------------------------------------------------------ phenoData <- pData(RGSet) phenoData[,1:6] ## ------------------------------------------------------------------------ manifest <- getManifest(RGSet) manifest head(getProbeInfo(manifest)) ## ------------------------------------------------------------------------ MSet <- preprocessRaw(RGSet) MSet ## ------------------------------------------------------------------------ head(getMeth(MSet)[,1:3]) head(getUnmeth(MSet)[,1:3]) ## ------------------------------------------------------------------------ RSet <- ratioConvert(MSet, what = "both", keepCN = TRUE) RSet ## ------------------------------------------------------------------------ beta <- getBeta(RSet) ## ------------------------------------------------------------------------ GRset <- mapToGenome(RSet) GRset ## ------------------------------------------------------------------------ beta <- getBeta(GRset) M <- getM(GRset) CN <- getCN(GRset) ## ------------------------------------------------------------------------ sampleNames <- sampleNames(GRset) probeNames <- featureNames(GRset) pheno <- pData(GRset) ## ------------------------------------------------------------------------ gr <- granges(GRset) head(gr, n = 3) ## ------------------------------------------------------------------------ annotation <- getAnnotation(GRset) names(annotation) ## ------------------------------------------------------------------------ head(getMeth(MSet)[,1:3]) head(getUnmeth(MSet)[,1:3]) ## ------------------------------------------------------------------------ qc <- getQC(MSet) qc plotQC(qc) ## ------------------------------------------------------------------------ densityPlot(MSet, sampGroups = phenoData$Sample_Group) ## ------------------------------------------------------------------------ densityBeanPlot(MSet, sampGroups = phenoData$Sample_Group) ## ------------------------------------------------------------------------ controlStripPlot(RGSet, controls="BISULFITE CONVERSION II") ## ---- eval=FALSE--------------------------------------------------------- # qcReport(RGSet, # sampNames = phenoData(RGSet)$Sample_Name, # sampGroups = phenoData(RGSet)$Sample_Group, # pdf= "qcReport.pdf") ## ------------------------------------------------------------------------ predictedSex <- getSex(GRset, cutoff = -2)$predictedSex predictedSex ## ------------------------------------------------------------------------ plotSex(getSex(GRset, cutoff = -2)) ## ------------------------------------------------------------------------ MSet.illumina <- preprocessIllumina(RGSet, bg.correct = TRUE, normalize = "controls") ## ------------------------------------------------------------------------ set.seed(666) MSet.swan <- preprocessSWAN(RGSet) ## ------------------------------------------------------------------------ GRset.quantile <- preprocessQuantile(RGSet, fixOutliers = TRUE, removeBadSamples = TRUE, badSampleCutoff = 10.5, quantileNormalize = TRUE, stratified = TRUE, mergeManifest = FALSE, sex = NULL) ## ------------------------------------------------------------------------ MSet.noob <- preprocessNoob(RGSet) ## ------------------------------------------------------------------------ GRset.funnorm <- preprocessFunnorm(RGSet) ## ---- eval = TRUE, echo = TRUE, mdsPlot---------------------------------- par(mfrow = c(2, 2)) # Following preprocessRaw() mdsPlot(getM(MSet.illumina), numPositions = 1000, sampGroups = targets$status, sampNames = targets$person, legendPos = "topright", main = "preprocessRaw") # Following preprocessQuantile() mdsPlot(getM(GRset.quantile), numPositions = 1000, sampGroups = targets$status, sampNames = targets$person, legendPos = "topright", main = "preprocessQuantile") # Following preprocessFunnorm() mdsPlot(getM(GRset.funnorm), numPositions = 1000, sampGroups = targets$status, sampNames = targets$person, legendPos = "bottomleft", main = "preprocessFunnorm") ## ------------------------------------------------------------------------ snps <- getSnpInfo(GRset) head(snps,10) ## ------------------------------------------------------------------------ GRset <- addSnpInfo(GRset) ## ------------------------------------------------------------------------ GRset <- dropLociWithSnps(GRset, snps=c("SBE","CpG"), maf=0) ## ---- eval=FALSE--------------------------------------------------------- # library(FlowSorted.Blood.450k) # cellCounts <- estimateCellCounts(RGSet) ## ------------------------------------------------------------------------ beta <- getBeta(GRset.funnorm) age <- pData(GRset.funnorm)$age dmp <- dmpFinder(beta, pheno = age, type = "continuous") head(dmp) ## ---- eval = TRUE, echo = TRUE, plotCpg---------------------------------- cpgs <- rownames(dmp)[1:4] par(mfrow = c(2, 2)) plotCpg(getBeta(GRset.funnorm), cpg = cpgs, pheno = age, type = "continuous", xlab = "age") ## ------------------------------------------------------------------------ pheno <- pData(GRset.funnorm)$status designMatrix <- model.matrix(~ pheno) designMatrix ## ---- eval=TRUE---------------------------------------------------------- bh_dmrs <- bumphunter(GRset.funnorm, design = designMatrix, cutoff = 0.2, B = 0, type = "Beta") ## ---- eval=FALSE--------------------------------------------------------- # # NOTE: Don't run, this takes a **long** time # bh_dmrs <- bumphunter(GRset.funnorm, # design = designMatrix, # cutoff = 0.2, # B = 1000, # type = "Beta") ## ---- eval=FALSE--------------------------------------------------------- # library(doParallel) # registerDoParallel(cores = 3) ## ---- eval=TRUE--------------------------------------------------------- names(bh_dmrs) head(bh_dmrs$table, n = 3) ## ---- eval = FALSE, echo = FALSE----------------------------------------- # # Ugly hack to make data available in vignette # # TODO: What's the proper way to do this? # lazyLoad(paste0(system.file(package = "tutorial.450k"), "/data/Rdata")) ## ---- eval=FALSE--------------------------------------------------------- # # NOTE: The `dmrs` object is included with the tutorial.450k package # library(tutorial.450k) # dmrs # head(dmrs$table) ## ---- eval = TRUE, echo = TRUE, Gviz------------------------------------- library(Gviz) genome <- "hg19" # NOTE: Using the non-permuted dmr <- bh_dmrs$table[1, ] chrom <- dmr$chr start <- dmr$start end <- dmr$end minbase <- start - 0.25 * (end - start) maxbase <- end + 0.25 * (end - start) pal <- c("#E41A1C", "#377EB8") # Start building the tracks iTrack <- IdeogramTrack(genome = genome, chromosome = dmr$chr, name = "") gTrack <- GenomeAxisTrack(col = "black", cex = 1, name = "", fontcolor = "black") # NOTE: This track takes a little while to create rTrack <- UcscTrack(genome = genome, chromosome = chrom, track = "refGene", from = minbase, to = maxbase, trackType = "GeneRegionTrack", rstarts = "exonStarts", rends = "exonEnds", gene = "name", symbol = "name2", transcript = "name", strand = "strand", fill = "darkblue", stacking = "squish", name = "RefSeq", showId = TRUE, geneSymbol = TRUE) # methylation data track gr <- granges(GRset.funnorm) gr$beta <- getBeta(GRset.funnorm) methTrack <- DataTrack(range = gr, groups = targets$status, genome = genome, chromosome = chrom, ylim = c(-0.05, 1.05), col = pal, type = c("a","p"), name = "DNA Meth.\n(beta value)", background.panel = "white", legend = TRUE, cex.title = 0.8, cex.axis = 0.8, cex.legend = 0.8) # DMR position data track dmrTrack <- AnnotationTrack(start = start, end = end, genome = genome, name = "DMR", chromosom = chrom) # Finally, plot the tracks tracks <- list(iTrack, gTrack, methTrack, dmrTrack, rTrack) sizes <- c(2, 2, 5, 2, 3) # set up the relative sizes of the tracks plotTracks(tracks, from = minbase, to = maxbase, showTitle = TRUE, add53 = TRUE, add35 = TRUE, grid = TRUE, lty.grid = 3, sizes = sizes, length(tracks)) ## ------------------------------------------------------------------------ library(sva) mval <- getM(GRset)[1:5000,] pheno <- pData(GRset) mod <- model.matrix(~as.factor(status), data=pheno) mod0 <- model.matrix(~1, data=pheno) sva.results <- sva(mval, mod, mod0) ## ---- eval=FALSE--------------------------------------------------------- # ab <- compartments(GRset.quantile, chr="chr14", resolution=100000) ## ------------------------------------------------------------------------ snps <- getSnpBeta(RGSet) head(snps) ## ------------------------------------------------------------------------ oob <- getOOB(RGSet) ## ---- eval = FALSE------------------------------------------------------- # library(shinyMethyl) # library(shinyMethylData) # runShinyMethyl(summary.tcga.raw, summary.tcga.norm) ## ------------------------------------------------------------------------ sessionInfo()