## ----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()