%\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{ELMER: Inferring Regulatory Element Landscapes and Transcription Factor Networks Using Methylomes} % !Rnw weave = knitr \documentclass{article} <>= BiocStyle::latex() @ \usepackage[utf8]{inputenc} \usepackage{graphicx} \usepackage{grffile} \usepackage{float} \usepackage{amsmath,epsfig,fullpage} \usepackage{hyperref} \usepackage{url} \usepackage[authoryear,round]{natbib} %\usepackage[OT1]{fonitenc} \usepackage{textcomp} %\\newcommand{\Rpackage}[1]{{\textit{#1}}} \title{ELMER: An R/Bioconductor Tool \\ Inferring Regulatory Element Landscapes and Transcription Factor Networks Using Methylomes} \author {Lijing Yao [cre, aut], Ben Berman [aut], Peggy Farnham [aut]\\ Hui Shen [ctb], Peter Laird [ctb], Simon Coetzee [ctb]} \begin{document} %\SweaveOpts{concordance=TRUE} % #' <>= % #' library(knitr) % #' opts_chunk$set(fig.align='center', par=TRUE) % #' knit_hooks$set(par=function(before, options, envir){ % #' if (before && options$fig.show!='none') par(mar=c(1,1,.1,.1),cex.lab=1,cex.axis=1,mgp=c(2,.7,0),tcl=-.3) % #' }, crop=hook_pdfcrop) % #' @ \maketitle \tableofcontents \section{Introduction} This document provides an introduction of the \Rpackage{ELMER}, which is designed to combine DNA methylation and gene expression data from human tissues to infer multi-level cis-regulatory networks. \Rpackage{ELMER} uses DNA methylation to identify enhancers, and correlates enhancer state with expression of nearby genes to identify one or more transcriptional targets. Transcription factor (TF) binding site analysis of enhancers is coupled with expression analysis of all TFs to infer upstream regulators. This package can be easily applied to TCGA public available cancer data sets and custom DNA methylation and gene expression data sets. ELMER analyses have 5 main steps: 1. Identify distal enhancer probes on HM450K.\\ 2. Identify distal enhancer probes with significantly different DNA methyaltion level in control group and experiment group.\\ 3. Identify putative target genes for differentially methylated distal enhancer probes.\\ 4. Identify enriched motifs for the distal enhancer probes which are significantly differentially methylated and linked to putative target gene.\\ 5. Identify regulatory TFs whose expression associate with DNA methylation at motifs.\\ \subsection{Installing and loading ELMER} To install this package, start R and enter <>= install.packages(devtools) library(devtools); devtools::install_github("lijingya/ELMER"); @ <>= source("http://bioconductor.org/biocLite.R") biocLite("ELMER") @ \section{Download example data} The following steps can be used to download the example data set for \Rpackage{ELMER}. <>= 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 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) @ \section{Quick start: running TCGA example} TCGA.pipe is a function for easily downloading TCGA data and performing all the analyses in ELMER. For illustration purpose, we skip the downloading step. The user can use the getTCGA function to download TCGA data or use TCGA.pipe by including "download" in the analysis option. The following command will do distal enhancer DNA methylation analysis and predict putative target genes, motif analysis and idnetify regulatory transcription factors. <>= 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"))) @ \section{Input data} MEE.data object is the input for all main functions of \Rpackage{ELMER}. MEE.data object can be generated by fetch.mee function. In order to perform all analyses in \Rpackage{ELMER}. MEE.data needs at least 4 files: a matrix of DNA methylation from HM450K platform for multiple samples; a matrix of gene expression for the same samples; a GRanges object containing the information for probes on HM450K such as names and coordinates; a gene annotation which is also a GRanges object. When TCGA data are used, the sample information will be automatically generated by fetch.mee function. However sample information should be provided when using custom data. \subsection{DNA methylation data} DNA methyaltion data feeding to \Rpackage{ELMER} should be a matrix of DNA methylation beta ($\beta$) value for samples (column) and probes (row) processed from row HM450K array data. If TCGA data were used, level 3 processed data from TCGA website will be downloaded and automatically transformed to the matrix by \Rpackage{ELMER}. The level 3 TCGA DNA methylation data were calculated as (M/(M+U)), where M represents the methylated allele intensity and U the unmethylated allele intensity. Beta values range from 0 to 1, reflecting the fraction of methylated alleles at each CpG in the each tumor; beta values close to 0 indicates low levels of DNA methylation and beta values close to 1 indicates high levels of DNA methylation. If user have raw HM450K data, these data can be processed by \Rpackage{Methylumi} or \Rpackage{minfi} generating DNA methylation beta ($\beta$) value for each CpG site and multiple samples. The getBeta function in \Rpackage{minfi} can be used to generate a matrix of DNA methylation beta ($\beta$) value to feed in \Rpackage{ELMER}. And we recommend to save this matrix as meth.rda since fetch.mee can read in files by specifying files' path which will help to reduce memory usage. <>= load("./ELMER.example/Result/LUSC/LUSC_meth_refined.rda") Meth[1:10, 1:2] @ \subsection{Gene expression data} Gene expresion data feeding to \Rpackage{ELMER} should be a matrix of gene expression values for samples (column) and genes (row). Gene expression value can be generated from different platforms: array or RNA-seq. The row data should be processed by other software to get gene or transcript level gene expression calls such as mapping by tophat, calling expression value by cufflink, RSEM or GenomeStudio for expression array. It is recommended to normalize expression data making gene expression comparable across samples such as quantile normalization. User can refer TCGA RNA-seq analysis pipeline to do generate comparable gene expression data. Then transform the gene expression values from each sample to the matrix for feeding into \Rpackage{ELMER}. If users want to use TCGA data, \Rpackage{ELMER} has functions to download the level 3 RNA-seq Hiseq data from TCGA website and transform the data to the matrix for feeding into \Rpackage{ELMER}. It is recommended to save this matrix as RNA.rda since fetch.mee can read in files by specifying the path of files which will help to reduce memory usage. <>= load("./ELMER.example/Result/LUSC/LUSC_RNA_refined.rda") GeneExp[1:10, 1:2] @ \subsection{Sample information} Sample information should be stored as a data.frame object containing sample ID, group labels (control and experiment). Sample ID and groups labels are required. Other information for each sample can be added to this data.frame object. When TCGA data were used, group labels (control and experiment) will be automatically generated by fetch.mee function by specifying option TCGA=TRUE. <>= mee <- fetch.mee(meth=Meth, exp=GeneExp, TCGA=TRUE) head(getSample(mee)) @ \subsection{Probe information} Probe information should be stored as a GRanges object containing the coordinates of each probe on the DNA methylation array and names of each probe. The default probe information is for HM450K fetching from \Rpackage{IlluminaHumanMethylation450kanno.ilmn12.hg19} <>= 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) @ \subsection{Gene information} Gene information should be stored as a GRanges object containing coordinates of each gene, gene id, gene symbol and gene isoform id. The default gene information is the UCSC gene annotation fetching from \Rpackage{Homo.sapiens} by \Rpackage{ELMER} function: txs. <>= 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) @ \subsection{MEE.data} MEE.data object is the input for multiple main functions of \Rpackage{ELMER}. MEE.data contains the above 5 components and making MEE.data object by fetch.mee function will keep each component consistent with each other. For example, althougth DNA methylation and gene expression matrixes have different rows (probe for DNA methylation and geneid for gene expression), the column (samples) order should be same in the two matrixes. fetch.mee function will keep them consistent when it generates the MEE.data object. <>= mee <- fetch.mee(meth=Meth, exp=GeneExp, TCGA=TRUE, probeInfo=probe, geneInfo=geneInfo) mee @ \section{Illustration of ELMER analysis steps} The example data set is a subset of chromosome 1 data from TCGA LUSC. ELMER analysis have 5 main steps which are shown below individually. TCGA.pipe introduced above is a pipeline combining all 5 steps and producing all results and figures. \subsection{Selection of probes within biofeatures} This step is to select HM450K probes, which locate far from TSS (at least 2Kb away) and within distal enhancer regions. These probes are called distal enhancer probes. For distal enhancer regions, we constructed a comprehensive list of putative enhancers combining REMC, ENCODE and FANTOM5 data. Enhancers were identified using ChromHMM for 98 tissues or cell lines from REMC and ENCODE project and we used the union of genomic elements labeled as EnhG1, EnhG2, EnhA1 or EnhA2 (representing intergenic and intragenic active enhancers) in any of the 98 cell types, resulting in a total of 389,967 non-overlapping enhancer regions. Additionally, FANTOM5 enhancers (43,011) identified by eRNAs for 400 distinct cell types were added to this list. Be default, this comprehensive list of putative enhancer and TSS annotated by GENCODE V15 and UCSC-gene will be used to select distal enhancer probes. But user can use their own TSS annotation or features such as H3K27ac ChIP-seq in a certain cell line. <>= #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") @ \subsection{Identifying differentially methylated probes} This step is to identify DNA methylation changes at distal enhancer probes which is carried out by function get.diff.meth. For each enhancer probe, function first ranked experiment samples and control samples by their DNA methylation beta values. To identify hypomethylated probes, function compared the lower control quintile (20\% of control samples with the lowest methylation) to the lower experiment quintile (20\% of experiment samples with the lowest methylation), using an unpaired one-tailed t-test. Only the lower quintiles were used because we did not expect all cases to be from a single molecular subtype, and we sought to identify methylation changes within cases from the same molecular subtype. 20\% (i.e. a quintile) was picked as a cutoff to include high enough sample numbers to yield t-test p-values that could overcome multiple hypothesis correction, yet low enough to be able to capture changes in individual molecular subtypes occurring in 20\% or more of the cases. This number can be set arbitrarily as an input to the get.diff.meth function in the \Rpackage{ELMER}, and should be tuned based on sample sizes in individual studies. The one tailed t-test was used to rule out the null hypothesis: $\mu experiment \geq \mu control$, where $\mu$experiment is the mean methylation within the lowest experiment quintile and $\mu$control is the mean within the lowest control quintile. Raw p-values were adjusted for multiple hypothesis testing using the Benjamini-Hochberg method, and probes were selected when they had adjusted p-value less than 0.01. For additional stringency, probes were only selected if the methylation difference: $\Delta = \mu experiment - \mu control$ was greater than 0.3. The same method was used to identify hypermethylated probes, except we used upper experiment quintile and upper control quintile, and chose the opposite tail in the t-test. If save parameter of get.diff.meth is true, two cvs files will be saved. If false, a data frame with the same content as the second file mentioned below will be reported. The first file contains all statistic results for all probes which were fed into the function. Based on this file, user can change different P value or sig.dir cutoff to select the significant results without redo the analysis. The second file contains statistic results for the probes that pass the significant criteria (P value and sig.dir). Both files contain four columns: probe, pvalue, ExperimentMinControl, adjust.p. 1. probe: the name of probes.\\ 2. pvalue: the raw P value from t-test.\\ 3. ExperimentMinControl: methylation difference $\Delta$.\\ 4. adjust.p: adjusted P value for t-test.\\ <>= ## 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") @ \subsection{Identifying putative probe-gene pairs} This step is to link enhancer probes with mehtylation changes to target genes with expression changes and report the putative target gene for selected probes. This is carried out by function get.pair. For each enhancer probe with differential methylation, the closest 10 upstream genes and the closest 10 downstream genes were tested for correlation between methylation of the probe and expression of the gene. To select these genes, the probe-gene distance was defined as the distance from the probe to a transcription start site specified by the UCSC gene annotation. Thus, exactly 20 statistical tests were performed for each probe, as follows. For each probe-gene pair, the samples (all experiment samples and control samples) were divided into two groups: the M group, which consisted of the upper methylation quintile (the 20\% of samples with the highest methylation at the enhancer probe), and the U group, which consisted of the lowest methylation quintile (the 20\% of samples with the lowest methylation.) The 20\%ile cutoff is a configurable parameter in the get.pair Default is 20\% as a balance, which would allow us to identify changes in a molecular subtype making up a minority (i.e. 20\%) of cases, while also yielding enough statistical power to make strong predictions. For each candidate probe-gene pair, the Mann-Whitney U test was used to test the null hypothesis that overall gene expression in group M was greater or equal than that in group U. This non-parametric test was used in order to minimize the effects of expression outliers, which can occur across a very wide dynamic range. For each probe-gene pair tested, the raw p-value Pr was corrected for multiple hypothesis using a permutation approach as follows (implemented in the get.permu function of the \Rpackage{ELMER} package). The gene in the pair was held constant, and x random methylation probes were used to perform the same one-tailed U test, generating a set of x permutation p-values (Pp). We chose the x random probes only from among those that were “distal” (greater than 2kb from an annotated transcription start site), in order to make these null-model probes qualitatively similar to the probe being tested. An empirical p-value Pe value was calculated using the following formula (which introduces a pseudo-count of 1): $$Pe = \frac{num(Pp \leq Pr)+ 1}{x+1}$$ This step is the most time consuming step since it requires a large amount calculations for permutation. The greater the permutation time is, the longer it will take. It is recommended to use multiple cores for this step. Default permutation time is 1000 which may need 12 hrs by 4 cores. However 10,000 permutations is recommended if high confidence results are desired but it may cost 2 days. If save paramter of get.pair function is true, two cvs files will be output. If save parameter is false, a data frame with the same contect as the second file mentioned as below will be output. The first file contains all statistic results for all probe-gene pairs. Based on this file, user can change different P value or sig.dir cutoff to select the significant results without redo the analysis. The second file contains statistic results for the probes that pass the significant criteria (Pe). Both files contain four columns: probe, GeneID, Symbol, Distance, Sides, Raw.p, Pe. 1. Probe: the name of probes.\\ 2. GeneID and Symbol is for the genes which are linked to the probe.\\ 3. Distance: the distance between the probe and the gene.\\ 4. Sides: right (R) side or left (L) side of probe and the rank based on distance. For example, L3 means the gene is the number 3 closest gene from the left side of the probe.\\ 5. Raw.p: P value from the Mann-Whitney U test for each pair.\\ 6. Pe: the empirical P value for each pair.\\ <>= ### 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") @ \subsection{Motif enrichment analysis on the selected probes} This step is to identify enriched motif in a set of probes which is carried out by function get.enriched.motif. The build in data Probes.motif is generated using FIMO with a p-value < 1e–4 to scan a +/- 100bp region around each probe using Factorbook motif position weight matrices (PWMs) and Jasper core human motif PWMs generated from the R package MotifDb. For each probe set tested (i.e. the list of gene-linked hypomethylated probes in a given experiment group), a motif enrichment Odds Ratio and a 95\% confidence interval were calculated using following formulas: $$ p= \frac{a}{a+b} $$ $$ P= \frac{c}{c+d} $$ $$ Odds\quad Ratio = \frac{\frac{p}{1-p}}{\frac{P}{1-P}} $$ $$ SD = \sqrt{\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d}} $$ $$ lower\quad boundary\quad of\quad 95\%\quad confidence\quad interval = \exp{(\ln{OR}-SD)} $$ where a is the number of probes within the selected probe set that contain one or more motif occurrences; b is the number of probes within the selected probe set that do not contain a motif occurrence; c and d are the same counts within the entire enhancer probe set. A probe set was considered significantly enriched for a particular motif if the 95\% confidence interval of the Odds Ratio was greater than 1.1 (specified by option lower.OR, 1.1 is default), and the motif occurred at least 10 times (specified by option min.incidence. 10 is default) in the probe set. As described in the text, Odds Ratios were also used for ranking candidate motifs. There will be two results if save parameter of get.enriched.motif is true. When save is false, only second result mentioned below will be reported. The first one is a csv file. This file contains the Odds Ratio and 95\% confidence interval for these Odds Ratios which pass the signficant cutoff (lower.OR and min.incidence). It contains 5 columns: motif, NumOfProbes, OR, lowerOR and upperOR. 1. motif: the name of motif.\\ 2. NumOfProbes: the number of probes with this motif in the given set of probes corresponding to min.incidence option.\\ 3. OR: the Odds Ratio.\\ 4. lowerOR: the lower boundary of 95\% confidence interval.\\ 5. upperOR: the upper boundary of 95\% confidence interval. The second file is a rda file listing each enriched motif and the probes containing the motif. <>= ### 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") @ \subsection{Identifying regulatory TFs} This step is to identify regulatory TF whose expression associates with TF binding motif DNA methylation which is carried out by function get.TFs.\\ For each motif considered to be enriched within a particular probe set, function will compare the average DNA methylation at all distal enhancer probes within +/- 100bp of a motif occurrence, to the expression of 1,982 human TFs A statistical test was performed for each motif-TF pair, as follows. The samples (all control and experiment samples) were divided into two groups: the M group, which consisted of the 20\% of samples with the highest average methylation at all motif-adjacent probes, and the U group, which consisted of the 20\% of samples with the lowest methylation. The 20th percentile cutoff is a parameter to the get.TFs function and was set to allow for identification of molecular subtypes present in 20\% of cases. For each candidate motif-TF pair, the Mann-Whitney U test was used to test the null hypothesis that overall gene expression in group M was greater or equal than that in group U. This non-parametric test was used in order to minimize the effects of expression outliers, which can occur across a very wide dynamic range. For each motif tested, this resulted in a raw p-value (Pr) for each of the 1982 TFs. All TFs were ranked by the -log10(Pr), and those falling within the top 5\% of this ranking were considered candidate upstream regulators. The best upstream TFs which are known recognized the motif was automatically extracted as putative regulatory TFs. If save parameter of get.TFs function is true, two files will be generated. If save parameter is false, only a data frame containing the same content with the first file mentioned below will be output. The first one csv file. This file contain the regulatory TF significantly associate with average DNA methylation at particular motif sites. It contains 4 columns: motif, top.potential.TF, potential.TFs and top\textunderscore 5percent. 1. motif: the names of motif.\\ 2. top.potential.TF: the highest ranking upstream TFs which are known recognized the motif.\\ 3. potential.TFs: TFs which are within top 5\% list and are known recognized the motif.\\ 4. top\textunderscore 5percent: all TFs which are within top 5\% list.\\ The second file is rda file. This file contains a matrix storing the statistic results for associations between TFs (row) and average DNA methylation at motifs (column). This matrix can be use to generate TF ranking plots by function TF.rank.plot. <>= ### 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") @ \section{Generating figures} \subsection{Scatter plots} Generate scatter plots for one probes' nearby 20 gene expression vs DNA methylation at this probe. Figure \ref{fig:figure1} <>= scatter.plot(mee,byProbe=list(probe=c("cg19403323"),geneNum=20), category="TN", save=TRUE) @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/cg19403323.byprobe.pdf} \protect\caption{SEach scatter plot shows the methylation level of an example probe cg19403323 in all LUSC samples plotted against the expression of one of 20 adjacent genes. \label{fig:figure1}} \end{center} \end{figure} \subsubsection{Scatter plot of one pair} Generate a scatter plot for one probe-gene pair. Figure \ref{fig:figure2} <>= scatter.plot(mee,byPair=list(probe=c("cg19403323"),gene=c("ID255928")), category="TN", save=TRUE,lm_line=TRUE) @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/cg19403323_SYT14.bypair.pdf} \protect\caption{Scatter plot shows the methylation level of an example probe cg19403323 in all LUSC samples plotted against the expression of the putative target gene SYT14. \label{fig:figure2}} \end{center} \end{figure} \subsubsection{TF expression vs. average DNA methylation} Generate scatter plot for TF expression vs average DNA methylation of the sites with certain motif. Figure \ref{fig:figure3} <>= 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) @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/TP53_TP63_TP73.byTF.pdf} \protect\caption{Each scatter plot shows the average methylation level of sites with the TP53 motif in all LUSC samples plotted against the expression of the transcription factor TP53, TP63, TP73 respectively. \label{fig:figure3}} \end{center} \end{figure} \subsection{Schematic plot} Schematic plot shows a breif view of linkages between genes and probes. To make a schematic plot, "Pair" object should be generated first. <>= # 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") @ \subsubsection{Nearby Genes} Generate schematic plot for one probe with 20 nearby genes and label the gene significantly linked with the probe in red. Figure \ref{fig:figure4} <>= schematic.plot(pair=pair, byProbe="cg19403323",save=TRUE) @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/cg19403323.schematic.byProbe.pdf} \protect\caption{The schematic plot shows probe colored in blue and the location of nearby 20 genes. The genes significantly linked to the probe were shown in red. \label{fig:figure4}} \end{center} \end{figure} \clearpage \subsubsection{Nearby Probes} Generate schematic plot for one gene with the probes which the gene is significantly linked to. Figure \ref{fig:figure5} <>= schematic.plot(pair=pair, byGene="ID255928",save=TRUE) @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/ID255928.schematic.byGene.pdf} \protect\caption{The schematic plot shows the gene colored in red and all blue colored probes, which are significantly linked to the expression of this gene. \label{fig:figure5}} \end{center} \end{figure} \subsection{Motif enrichment plot} Motif enrichment plot shows the enrichment levels for the selected motifs. Figure\ref{fig:figure6} <>= 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. @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/hypo.motif.enrichment.pdf} \protect\caption{The plot shows the Odds Ratio (x axis) for the selected motifs with OR above 1.3 and lower boundary of OR above 1.3. The range shows the 95\% confidence interval for each Odds Ratio. \label{fig:figure6}} \end{center} \end{figure} \clearpage \subsection{TF ranking plot} TF ranking plot shows statistic -log10(P value) assessing the anti-correlation level of TFs expression level with average DNA methylation level at sites with a given motif. Figure \ref{fig:figure7} <>= 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) @ \begin{figure}[H] \begin{center} \includegraphics[scale=0.5]{ELMER.example/TP53.TFrankPlot.pdf} \protect\caption{Shown are TF ranking plots based on the score (-log(P value)) of association between TF expression and DNA methylation of the TP53 motif in the LUSC cancer type . The dashed line indicates the boundary of the top 5\% association score. The top 3 associated TFs and the TF family members (dots in red) that are associated with that specific motif are labeled in the plot. \label{fig:figure7}} \end{center} \end{figure} \newpage <>= sessionInfo() @ \end{document}