## ----style, eval=TRUE, echo=FALSE, results="asis"------------------------ BiocStyle::latex() ## ----installing, eval=FALSE,hide=TRUE, message=FALSE,include=FALSE------- ## install.packages(devtools) ## library(devtools); ## devtools::install_github("lijingya/ELMER"); ## ----install, eval=FALSE------------------------------------------------- ## source("http://bioconductor.org/biocLite.R") ## biocLite("ELMER") ## ----example.data.run, echo=FALSE, hide=TRUE, message=FALSE, include=FALSE---- if(!file_test("-d", "./ELMER.example")) { if(Sys.info()["sysname"] == "Windows") mode <- "wb" else mode <- "w" URL <- "https://www.dropbox.com/s/mxf2u1wcauenx8m/ELMER.example.zip?raw=1" downloader::download(URL,destfile = "ELMER.example.zip",mode= mode) utils::unzip("ELMER.example.zip") } library(ELMER, quietly = TRUE) ## ----example.data, eval=FALSE-------------------------------------------- ## #Example file download from URL: https://dl.dropboxusercontent.com/u/61961845/ELMER.example.tar.gz ## URL <- "https://dl.dropboxusercontent.com/u/61961845/ELMER.example.tar.gz" ## download.file(URL,destfile = "ELMER.example.tar.gz",method= "curl") ## untar("./ELMER.example.tar.gz") ## library(ELMER) ## ----tcga.pipe,eval=FALSE------------------------------------------------ ## TCGA.pipe("LUSC",wd="./ELMER.example",cores=detectCores()/2,permu.size=300,Pe=0.01, ## analysis = c("distal.probes","diffMeth","pair","motif","TF.search"), ## diff.dir="hypo",rm.chr=paste0("chr",c("X","Y"))) ## ----dna.methylation.data------------------------------------------------ load("./ELMER.example/Result/LUSC/LUSC_meth_refined.rda") Meth[1:10, 1:2] ## ----gene.expression.data------------------------------------------------ load("./ELMER.example/Result/LUSC/LUSC_RNA_refined.rda") GeneExp[1:10, 1:2] ## ----sample.information-------------------------------------------------- mee <- fetch.mee(meth=Meth, exp=GeneExp, TCGA=TRUE) head(getSample(mee)) ## ----probe.information--------------------------------------------------- probe <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19, what="Locations") probe <- GRanges(seqnames=probe$chr, ranges=IRanges(probe$pos, width=1, names=rownames(probe)), strand=probe$strand, name=rownames(probe)) mee <- fetch.mee(meth=Meth, exp=GeneExp, TCGA=TRUE, probeInfo=probe) getProbeInfo(mee) ## ----gene.information---------------------------------------------------- geneAnnot <- txs() ## In TCGA expression data, geneIDs were used as the rowname for each row. However, numbers ## can't be the rownames, "ID" was added to each gene id functioning as the rowname. ## If your geneID is consistent with the rownames of the gene expression matrix, adding "ID" ## to each geneID can be skipped. geneAnnot$GENEID <- paste0("ID",geneAnnot$GENEID) geneInfo <- promoters(geneAnnot,upstream = 0, downstream = 0) save(geneInfo,file="./ELMER.example/Result/LUSC/geneAnnot.rda") mee <- fetch.mee(meth=Meth, exp=GeneExp, TCGA=TRUE, geneInfo=geneInfo) getGeneInfo(mee) ## ----mee.data------------------------------------------------------------ mee <- fetch.mee(meth=Meth, exp=GeneExp, TCGA=TRUE, probeInfo=probe, geneInfo=geneInfo) mee ## ----selection.of.probes.within.biofeatures------------------------------ #get distal enhancer probes that are 2kb away from TSS and overlap with REMC and FANTOM5 #enhancers on chromosome 1 Probe <- get.feature.probe(rm.chr=paste0("chr",c(2:22,"X","Y"))) save(Probe,file="./ELMER.example/Result/LUSC/probeInfo_feature_distal.rda") ## ----identifying.differentially.methylated.probes, eval=FALSE------------ ## ## fetch.mee can take path as input. ## mee <- fetch.mee(meth="./ELMER.example/Result/LUSC/LUSC_meth_refined.rda", ## exp="./ELMER.example/Result/LUSC/LUSC_RNA_refined.rda", TCGA=TRUE, ## probeInfo="./ELMER.example/Result/LUSC/probeInfo_feature_distal.rda", ## geneInfo="./ELMER.example/Result/LUSC/geneAnnot.rda") ## ## sig.diff <- get.diff.meth(mee, cores=detectCores()/2, dir.out ="./ELMER.example/Result/LUSC", ## diff.dir="hypo", pvalue = 0.01) ## ## ## sig.diff[1:10,] ## significantly hypomethylated probes ## ## # get.diff.meth automatically save output files. ## # getMethdiff.hypo.probes.csv contains statistics for all the probes. ## # getMethdiff.hypo.probes.significant.csv contains only the significant probes which ## # is the same with sig.diff ## dir(path = "./ELMER.example/Result/LUSC", pattern = "getMethdiff") ## ----identifying.putative.probe.gene.pairs, eval=FALSE------------------- ## ### identify target gene for significantly hypomethylated probes. ## ## Sig.probes <- read.csv("./ELMER.example/Result/LUSC/getMethdiff.hypo.probes.significant.csv", ## stringsAsFactors=FALSE)[,1] ## head(Sig.probes) # significantly hypomethylated probes ## ## ## Collect nearby 20 genes for Sig.probes ## nearGenes <-GetNearGenes(TRange=getProbeInfo(mee,probe=Sig.probes), ## geneAnnot=getGeneInfo(mee),cores=detectCores()/2) ## ## ## Identify significant probe-gene pairs ## Hypo.pair <-get.pair(mee=mee,probes=Sig.probes,nearGenes=nearGenes, ## permu.dir="./ELMER.example/Result/LUSC/permu",permu.size=200,Pe = 0.01, ## dir.out="./ELMER.example/Result/LUSC",cores=detectCores()/2,label= "hypo") ## ## head(Hypo.pair) ## significant probe-gene pairs ## ## # get.pair automatically save output files. ## #getPair.hypo.all.pairs.statistic.csv contains statistics for all the probe-gene pairs. ## #getPair.hypo.pairs.significant.csv contains only the significant probes which is ## # same with Hypo.pair. ## dir(path = "./ELMER.example/Result/LUSC", pattern = "getPair") ## ----motif.enrichment.analysis.on.selected.probes, eval=FALSE------------ ## ### identify enriched motif for significantly hypomethylated probes which ## ##have putative target genes. ## ## Sig.probes.paired <- read.csv("./ELMER.example/Result/LUSC/getPair.hypo.pairs.significant.csv", ## stringsAsFactors=FALSE)[,1] ## head(Sig.probes.paired) # significantly hypomethylated probes with putative target genes ## ## enriched.motif <-get.enriched.motif(probes=Sig.probes.paired, ## dir.out="./ELMER.example/Result/LUSC", label="hypo", ## min.incidence = 10,lower.OR = 1.1) ## names(enriched.motif) # enriched motifs ## head(enriched.motif["TP53"]) ## probes in the given set that have TP53 motif. ## ## # get.enriched.motif automatically save output files. ## # getMotif.hypo.enriched.motifs.rda contains enriched motifs and the probes with the motif. ## # getMotif.hypo.motif.enrichment.csv contains summary of enriched motifs. ## dir(path = "./ELMER.example/Result/LUSC", pattern = "getMotif") ## ## # motif enrichment figure will be automatically generated. ## dir(path = "./ELMER.example/Result/LUSC", pattern = "motif.enrichment.pdf") ## ----identifying.regulatory.tf, eval=FALSE------------------------------- ## ### identify regulatory TF for the enriched motifs ## ## load("./ELMER.example/Result/LUSC/getMotif.hypo.enriched.motifs.rda") ## TF <- get.TFs(mee=mee, enriched.motif=enriched.motif,dir.out="./ELMER.example/Result/LUSC", ## cores=detectCores()/2, label= "hypo") ## ## # get.TFs automatically save output files. ## # getTF.hypo.TFs.with.motif.pvalue.rda contains statistics for all TF with average ## # DNA methylation at sites with the enriched motif. ## # getTF.hypo.significant.TFs.with.motif.summary.csv contains only the significant probes. ## dir(path = "./ELMER.example/Result/LUSC", pattern = "getTF") ## ## # TF ranking plot based on statistics will be automatically generated. ## dir(path = "./ELMER.example/Result/LUSC/TFrankPlot", pattern = "pdf") ## ----eval=FALSE---------------------------------------------------------- ## scatter.plot(mee,byProbe=list(probe=c("cg19403323"),geneNum=20), ## category="TN", save=TRUE) ## ----eval=FALSE---------------------------------------------------------- ## scatter.plot(mee,byPair=list(probe=c("cg19403323"),gene=c("ID255928")), ## category="TN", save=TRUE,lm_line=TRUE) ## ----eval=FALSE---------------------------------------------------------- ## load("ELMER.example/Result/LUSC/getMotif.hypo.enriched.motifs.rda") ## scatter.plot(mee,byTF=list(TF=c("TP53","TP63","TP73"), ## probe=enriched.motif[["TP53"]]), category="TN", ## save=TRUE,lm_line=TRUE) ## ----eval=FALSE---------------------------------------------------------- ## # Make a "Pair" object for schematic.plot ## pair <- fetch.pair(pair="./ELMER.example/Result/LUSC/getPair.hypo.pairs.significant.withmotif.csv", ## probeInfo = "./ELMER.example/Result/LUSC/probeInfo_feature_distal.rda", ## geneInfo = "./ELMER.example/Result/LUSC/geneAnnot.rda") ## ----eval=FALSE---------------------------------------------------------- ## schematic.plot(pair=pair, byProbe="cg19403323",save=TRUE) ## ----eval=FALSE---------------------------------------------------------- ## schematic.plot(pair=pair, byGene="ID255928",save=TRUE) ## ----eval=FALSE---------------------------------------------------------- ## motif.enrichment.plot(motif.enrichment="./ELMER.example/Result/LUSC/getMotif.hypo.motif.enrichment.csv", ## significant=list(OR=1.3,lowerOR=1.3), ## label="hypo", save=TRUE) ## different signficant cut off. ## ----eval=FALSE---------------------------------------------------------- ## load("./ELMER.example/Result/LUSC/getTF.hypo.TFs.with.motif.pvalue.rda") ## TF.rank.plot(motif.pvalue=TF.meth.cor, motif="TP53", TF.label=list(TP53=c("TP53","TP63","TP73")), ## save=TRUE) ## ----finalsession-------------------------------------------------------- sessionInfo()