--- title: "**Analyzing MeRIP-seq data with TRESS**" author: - name: Zhenxing Guo email: zhenxing.guo@emory.edu - name: Hao Wu email: hao.wu@emory.edu package: TRESS abstract: > This vignette introduces the use of the package TRESS (${\it T}$oobox for m${\it R}$NA ${\it E}$pigenetics ${\it S}$equencing analysi${\it S}$), which is designed for analysis of MeRIP-seq data. It provides functionalities for m6A methylation identification. TRESS utilizes Bayesian hierarchical negative binomial models to extract signals in raw read counts, and then conducts Wald tests to detect m6A methylation. vignette: > %\VignetteIndexEntry{The TRESS User's Guide} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} output: BiocStyle::html_document: toc_float: FALSE --- \tableofContents # Introduction ## Background The post-transcriptional epigenetic modification on mRNA is an emerging field to study the gene regulatory mechanism and their association with diseases. Recently developed high-throughput sequencing technology named Methylated RNA Immunoprecipitation Sequencing (MeRIP-seq) enables one to profile mRNA epigenetic modification transcriptome-wide. In MeRIP-seq, mRNA is first fragmented into approximately 100-nucleotide-long oligonucleotides, and then immunoprecipitated by an anti-m$^6$A affinity purified antibody. In addition to the immunoprecipitated (IP) samples, libraries are also prepared for input control fragments to measure the corresponding reference mRNA abundance. This process is an RNA-seq experiment. After sequencing, the reads from both the IP and the input samples are aligned to the reference genome. Due to the enrichment from IP process, transcriptomic regions with m$^6$A will have more reads clustered and have peak-like shapes when visualizing the read counts along the genome. Therefore, people often refer the m$^6$A regions as "peaks", which is a term usually used in ChIP-seq to represent the protein binding sites. Figure \@ref(fig:m6aPeak) shows some example peaks on the Fat2 gene from a dataset to study m$^6$A dynamics during mouse brain development, where m$^6$A in cerebellums from 2-week old mice are profiled with two biological replicates. ```{r m6aPeak, echo=FALSE, out.height = "80%", out.width = "80%", include=TRUE,fig.cap="Peaks from MeRIP-seq in a mouse brain study"} knitr::include_graphics("m6APeak.png") ``` Two major tasks in the analysis of MeRIP-seq data are to identify transcriptome-wide m$^6$A methylation (namely "peak calling"), and differential m6A methylation. For the first problem, TRESS builds a two-step procedure to identify transcriptome wide m$^6$A regions. In the first step, it quickly scans the whole transcriptome and losely identify candidate regions using an ad hoc algorithm. In the second step, it models read counts from candidate regions using an empirical hierarchical negative binomial model to accounts for all-sources variations. It also imposes a prior on the dispersion of methylation, which induces a shrinkage variance estimate by borrowing information from all candidate regions. Wald test is constructed to detect significant m$^6$A regions from the candidates. Metthod for the second problem is under development and will be available soon. ## Installation From GitHub: ```{r, eval = FALSE, warning=FALSE, message=FALSE} install.packages("devtools") # if you have not installed "devtools" package library(devtools) install_github("https://github.com/ZhenxingGuo0015/TRESS", build_vignettes = TRUE) ``` To view the package vignette in HTML format, run the following lines in R ```{r eval=FALSE,warning=FALSE,message=FALSE} library(TRESS) vignette("TRESS") ``` # Input data preparation{#section:InputData} TRESS requires paired input control and IP BAM files for each replicate of all samples: "input1.bam \& ip1.bam", "input2.bam \& ip2.bam", .... in order to call peaks from all replicates. The BAM files contain mapped reads sequenced from respective samples and are the output of sequence alignment tools like ``Bowtie2``. For illustration purpose, we include four example BAM files and one corresponding genome annotation file in our publicly available data package ``datasetTRES``on github, which can be installed with ```{r, eval= FALSE} install_github("https://github.com/ZhenxingGuo0015/datasetTRES") ``` The BAM files contain sequencing reads (only on chromosome 19) from two input \& IP mouse brain cerebellum samples. In addition to BAM files, TRESS also needs a path to an annotation file in order to obtain transcriptome-wide bins, bin-level read counts and annotation of each peak. The annotation file is actually is a TXDB and is saved in format of ``*.sqlite``, which is easily created using R function ``makeTxDbFromUCSC()`` from Bioconductor package ``GenomicFeatures``: ```{r, eval=FALSE, message=FALSE, warning=FALSE} ## Directly use "makeTxDbFromUCSC" function to create one library(GenomicFeatures) txdb = makeTxDbFromUCSC("mm9", "knownGene") # saveDb(txdb, file = paste0("YourPATH", "/", "YourGenome.sqlite") ## or load a TxDb annotation package like # library(TxDb.Mmusculus.UCSC.mm9.knownGene) # txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene # saveDb(txdb, file = paste0("YourPATH", "/", "mm9_knownGene.sqlite") ``` If one already has one .gtf file, then a TXDB file can be created using the following code: ```{r, eval=FALSE, message=FALSE, warning=FALSE} library(GenomicFeatures) txdb = makeTxDbFromGFF("directory/to/your/xxx.gtf", format = "gtf") # saveDb(txdb, file = paste0("YourPATH", "/", "YourGenome.sqlite") ``` # Example dataset {#section:EData} TRESS provides one example MeRIP-seq dataset. The raw sequencing reads were downloaded from GEO, and mapped using ``Bowtie2`` to generate corresponding BAM files. With BAM files, we apply TRESS to obtain transcriptome bins and candidate regions. A small subset of both bin-level and region-level data are stored as examples in TRESS for package illustration purpose. The example dataset ``Basal`` comes from 7 basal mouse brain samples (GEO accession number GSE113781). The total number of bins and candidate regions for original data are 1,620,977 and 13,017 respectively. Due to space limit, only data of randomly selected 1000 bins and 500 regions (bins and regions are not necessarily overlapped) are saved. In particular, for each bin and candidate region, both the genomic coordinate and read counts (across 7 paired input \& ip replicates) are kept. ```{r loadData,warning=FALSE,message=FALSE} library(TRESS) data("Basal") Bins <- Basal$Bins$Bins BinCounts <- Basal$Bins$Counts dim(BinCounts) Candidates <- Basal$Candidates$Regions CandidatesCounts <- Basal$Candidates$Counts sf <- Basal$Bins$sf ``` Check the coordinates of each bin and candidate region. ```{r BasalBins} head(Bins, 3) head(Candidates, 3) ``` Check the read counts in each bin and candidate region. ```{r BasalCount} dim(BinCounts) head(BinCounts, 3) dim(CandidatesCounts) head(CandidatesCounts,3) ``` Check the library size factor of each sample in this data. ```{r Basalsf} sf ``` # Detection of transcriptome wide peaks {#section:Peak} Peak calling in TRESS is performed using one wrapper function ``TRESS_peak()``. Here we use example BAM files from datasetTRES as an example. ```{r, eval= FALSE,warning=FALSE, message=FALSE} ## Directly take BAM files in "datasetTRES" available on github library(datasetTRES) Input.file = c("cb_input_rep1_chr19.bam", "cb_input_rep2_chr19.bam") IP.file = c("cb_ip_rep1_chr19.bam", "cb_ip_rep2_chr19.bam") BamDir = file.path(system.file(package = "datasetTRES"), "extdata/") annoDir = file.path(system.file(package = "datasetTRES"), "extdata/mm9_chr19_knownGene.sqlite") OutDir = "/directory/to/output" TRESS_peak(IP.file = IP.file, Input.file = Input.file, Path_To_AnnoSqlite = annoDir, InputDir = BamDir, OutputDir = OutDir, # specify a directory for output experiment_name = "examplebyBam", # name your output filetype = "bam") ``` ```{r, eval= TRUE} ### example peaks peaks = read.table(file.path(system.file(package = "TRESS"), "extdata/examplebyBam_peaks.xls"), sep = "\t", header = TRUE) head(peaks[, -c(5, 14, 15)], 3) ``` To replace the example BAM files with your BAM files, the codes are: ```{r, eval=FALSE} ## or, take BAM files from your path Input.file = c("input_rep1.bam", "input_rep2.bam") IP.file = c("ip_rep1.bam", "ip_rep2.bam") BamDir = "/directory/to/BAMfile" annoDir = "/path/to/xxx.sqlite" OutDir = "/directory/to/output" TRESS_peak(IP.file = IP.file, Input.file = Input.file, Path_To_AnnoSqlite = annoDir, InputDir = BamDir, OutputDir = OutDir, experiment_name = "xxx", filetype = "bam") peaks = read.table(paste0(OutDir, "/", "xxx_peaks.xls"), sep = "\t", header = TRUE) head(peaks, 3) ``` As a wrapper function, ``TRESS_peak()`` combines multiple subfunctions to detect m6A regions transcriptome-wide. The whole process involves following steps: $\quad 1$. Divide genome into equal sized bins and obtain read counts in each replicate using ``DivideBins()``. $\quad 2$. Call candidate regions based on bin-level data across all replicates using ``CallCandidates()`` $\quad 3$. Fit Negative binomial models for read counts in candidate regions to detect and rank peaks among all candidates using ``CallPeaks.multiRep()``. If there is only one replicate, functions in step 2 and step 3 will be replaced by one function ``CallPeaks.oneRep()``, while function ``DivideBins()`` will still be used to obtain bin-level data. Please see the man pages for detailed usage of ``DivideBins()``. Given bin-level data, peak calling by sub-functions are included in following Section \@ref(section:PeakMultiReps) and Section \@ref(section:PeakOneRep). ## Peak calling with multiple replicates {#section:PeakMultiReps} ### Obtain candidate regions With bin-level read counts, binomial test is first conducted for each bin. Then an ad hoc bump-finding algorithm is applied to merge significant bins (and/or bins with large fold change) and form bumps in each replicate. Bumps from all input \& IP replicates are unioned together to construct a list of candidate regions. Both binomial tests and bump-finding are done by function ``CallCandidates()`` in TRESS. Here, we use example dataset ``Basal`` introduced in Section \@ref(section:EData) to run ``CallCandidates()``. ```{r,eval=TRUE, message=FALSE, warning=FALSE} ## load in first example dataset data("Basal") Candidates = CallCandidates( Counts = Basal$Bins$Counts, bins = Basal$Bins$Bins ) ``` ### Call peaks from candidates After obtaining candidate regions, TRESS models read counts in candidate regions using a hierarchical negative-binomial model. Wald tests are then conducted to detect significant regions from candidates, where a region with methylation level significantly higher than the background is considered as sigificant. The background methylation level is estimated based on total read counts from non-candidate regions, which is calculated as the total bin-counts subtracted by counts from candidate regions. This can be obtained using function ``BgMethylevel()`` in TRESS. With background methylation level estimated, complete parameter estimation and statistical inference for candidate regions are achieved by function ``CallPeaks.multiRep()`` in TRESS, where argument ``Candidates`` is a list containing at least two components: genomic coordinates (e.g., "chr", "start" and "end" etc) of all candidate regions and their read counts across all samples. Here is an example, provided that the background methylation level is 0.5, ```{r, eval=TRUE, message= TRUE, warning= FALSE} data("Basal") ### load candidate regions Basal$Candidates$sf = Basal$Bins$sf peaks1 = CallPeaks.multiRep( Candidates = Basal$Candidates, mu.cutoff = 0.5 ) head(peaks1, 3) ``` _**A word about CallPeaks.multiRep()**_ Different criteria are available in TRESS to filter peaks based on results of Wald tests, including p-values, FDR and log fold change (logFC). One needs to first specify a criterion through argument ``WhichThreshold`` and then provide a cutoff for that criterion through one or two arguments named ``*.cutoff``. By default, TRESS adopts both FDR and logFC (through``WhichThreshold = "fdr_lfc"``) to select peaks. The default cutoffs for FDR and logFC are 0.05 and 0.7 (for fold change of 2) respectively, meaning that peaks with FDR < 0.05 and logFC > 0.7 will be kept. With ``WhichThreshold = "fdr_lfc"``, one can change the cutoffs to more stringent values, e.g., FDR < 0.01 and logFC > 1.6 by setting ``fdr.cutoff = 0.01`` and ``lfc.cutoff = 1.6``: ```{r, eval=TRUE, message= TRUE, warning= FALSE} ### use different threshold to filter peaks peaks2 = CallPeaks.multiRep( Candidates = Basal$Candidates, mu.cutoff = 0.5, fdr.cutoff = 0.01, lfc.cutoff = 1.6 ) ``` However, one can also set ``WhichThreshold = "fdr"`` to select peaks only with FDR, then by default peaks with FDR < 0.05 will be kept: ```{r, eval=TRUE, message= TRUE, warning= FALSE} ### use different threshold to filter peaks peaks3 = CallPeaks.multiRep( Candidates = Basal$Candidates, mu.cutoff = 0.5, WhichThreshold = "fdr" ) ``` Please read the manual page of function ``CallPeaks.multiRep()`` for more details on each argument. Given the usage of function "CallPeaks.multiRep()", it can also be adopted to re-rank existing peaks with our developed methods, if their read counts are available. This may perform bad if one didn't properly estimate library size factor (one argument of ``Candidates``) for each sample. Based on our experience, the estimation of size factor should be based on the bin-level counts across the whole transcriptome, not the region-level counts. For background methylation level, one can use 0.5 but it would be informative if estimating it from the data. ## Peak calling with one replicate {#section:PeakOneRep} The above two-step approach is for data with multiple replicates. For data with only one replicate, ``TRESS_peak()`` calls function ``CallPeaks.oneRep()`` for peak calling given bin-level data. In this case, bumps in the only one replicate are taken as final list of peaks. The statistical significance of each peak comes from binomial tests. Here is an example of how to run function ``CallPeaks.oneRep()``. ```{r, eval = TRUE, message= FALSE, warning= FALSE} # A toy example data("Basal") bincounts = Basal$Bins$Counts[, 1:2] sf0 = Basal$Bins$sf[1:2] bins = Basal$Bins$Bins peaks = CallPeaks.oneRep(Counts = bincounts, sf = sf0, bins = bins) head(peaks, 3) ``` ## Peak visualization With pre-called peaks in hand, one can visualize them using function "ShowOnePeak()" in TRESS. The usage of this function is ```{r, eval=FALSE, message= FALSE, warning= FALSE} ShowOnePeak(onePeak, allBins, binCounts, ext = 500, ylim = c(0,1)) ``` In order to run this function, one needs to have: 1) "onePeak": a pre-called peak saved as a dataframe, which contains genomic positions for that peak: "chr", "start", "end", "strand"; 2) "allBins": genomic positions ("chr", "start", "end", "strand") of bins overlapping at least with the pre-called peak; 3) "binCounts": the corresponding bin-level read counts in each replicate. This function will plot for each replicate: the methylation level of bins (blue bars) within the target peak (shade region in pink), and the normalized sequencing depth for input samples (curves in grey). We show some example plots here: ```{r, eval=TRUE, message= FALSE, warning= FALSE} peaks = read.table(file.path(system.file(package = "TRESS"), "extdata/examplebyBam_peaks.xls"), sep = "\t", header = TRUE) load(file.path(system.file(package = "TRESS"), "extdata/examplebyBam.rda")) allBins = as.data.frame(bins$bins) colnames(allBins)[1] = "chr" allBins$strand = binStrand head(peaks, 1) ShowOnePeak(onePeak = peaks[1,], allBins = allBins, binCounts = allCounts) ``` # Session info {.unnumbered} ```{r sessionInfo, echo=FALSE} sessionInfo() ```