%\VignetteIndexEntry{Identify differential APA usage from RNA-seq alignments} %\VignettePackage{roar} \documentclass{article} \author{Elena Grassi\\ Department of Molecular Biotechnologies and Health Sciences,\\ MBC, University of Turin, Italy\\ } \usepackage{Sweave} \usepackage{hyperref} % From dexseq vignette \newcommand{\Robject}[1]{{\small\texttt{#1}}} \newcommand{\Rfunction}[1]{\Robject{#1}} \newcommand{\Rpackage}[1]{\textit{#1}} \newcommand{\Rclass}[1]{\textit{#1}} \title{Identify differential APA usage from RNA-seq alignments} \date{\Rpackage{roar} version \Sexpr{packageDescription("roar")$Version} (Last revision 2014-09-24)} \begin{document} \SweaveOpts{concordance=TRUE} \maketitle \begin{abstract} This vignette describes how to use the Bioconductor package \Rpackage{roar} to detect preferential usage of shorter isoforms via alternative poly-adenylation from RNA-seq data. The approach is based on Fisher test to detect disequilibriums in the number of reads falling over the 3'UTRs when comparing two biological conditions. The name roar means ``ratio of a ratio'': counts and fragments lengths are used to calculate the prevalence of the short isoform over the long one in both conditions, therefore the ratio of these ratios represents the relative ``shortening'' (or lengthening) in one condition with respect to the other. As input, \Rpackage{roar} uses alignments files for the two conditions and coordinates of the 3'UTRs with alternative polyadenylation sites. Here, the method is demonstrated on the data from the package \Rpackage{RNAseqData.HNRNPC.bam.chr14}. To cite this software, please refer to \Rfunction{citation("roar")}. \end{abstract} \tableofcontents \section{Introduction} The alternative polyadenylation mechanism at the basis of the existence of short and long isoforms of the same gene is reviewed in Elkon et al~\cite{Elkon2013}. For the relevance of this phenomenon to development and cancer see~\cite{Sandberg2008,Mayr2009,Ji2009}. \section{Input data} A roar analysis starts from some alignments, obtained via standard RNAseq experiments and data processing. It is possible to build a \Rclass{RoarDataset} object giving its constructor two lists of bam files obtained from samples of the two conditions to be compared or to use another constructor that takes directly two lists of \Rclass{GAlignments}. \subsection{Alternative polyadenylation annotations} The other information needed to build a \Rclass{RoarDataset} object are 3'UTRs coordinates (with canonical and alternative polyadenylation sites) for the genes that one wants to analyze - these could be given using a gtf file or a GRanges object. The gtf file should have an attribute (\emph{metadata column} for the \Rclass{GRanges} object) called ``gene\_id'' that ends with ``\_PRE'' or ``\_POST'' to address respectively the short and the long isoforms. An element in the annotation is considered ``PRE'' (i.e. common to the short and long isoform of the transcript) if its gene\_id ends with ``\_PRE''. If it ends with ``\_POST'' it is considered the portion present only in the long isoform. The prefix of gene\_id should be a unique identifier for the gene and each identifier has to be associated with only one ``\_PRE'' and one ``\_POST'', leading to two genomic regions associated to each gene\_id. The gtf can also contain an attribute (\emph{metadata column} for the \Rclass{GRanges} object) that represents the lengths of PRE and POST portions on the transcriptome. If this is omitted the lengths on the genome are used instead (this will results in somewhat inaccurate length corrections for genes with either the PRE or the POST portion comprising an intron). Note that right now every gtf entry (or none of them) should have it. In the package we added a file (examples/apa.gtf) that follows these directives: it is build upon the hg19 human genome release using PolyA\_DB (\cite{Cheng2006}) version 2 to obtain coordinates of alternative polyadenylation (from now on addressed as APA) sites. Every gene with at least one APA site is considered; their longest transcript is used to get the ``classical'' polyadenylation syte (canonical end of the longest transcript stored in UCSC) and the most proximal (with respect to the TSS, that is the farther to the canonical end of the transcript; but preferring sites still in the 3'UTR) reported APA site found in PolyA\_DB is used to define the end of the short isoform. Therefore we define the ``PRE'' portion as the stretch of DNA starting from the beginning of the exon containing (or proximal to) the APA site and ending at the site position, while the POST portion starts there and ends at the canonical transcript end. Using only the nearest exon to the APA site and not the entire transcript to define the short isoform avoids spurious signals deriving from other alternative splicing events. To summarize the package requirements are: every entry in the gtf should have an attribute called ``gene\_id'' formed by a prefix representing unique gene identifiers and a suffix. The suffix defines if the given coordinates refer to the portion of this gene common between the short and the long isoform (``\_PRE'') or to the portion pertaining only to the long isoform (``\_POST''). Every gene identifier must appear two (and only two) times in the gtf, one with the suffix ``\_PRE'' and one with ``\_POST''. An example of how PRE and POST are defined is reported in Figure~\ref{Fig1}. \begin{figure}[h] \includegraphics[width=\textwidth]{prepost.pdf} \caption{An example of PRE/POST portions definition} \label{Fig1} \end{figure} \section{Analysis steps} The suggestion is to use one of the two wrapper scripts (\emph{roarWrapper} or \emph{roarWrapper\_chrBychr}) or at least follow their tracks. The principal steps of the analysis are performed by different methods that receive a \Rclass{RoarDataset} object as an argument and returns it with the needed step performed. It is also possible to call directly the methods to get the results (eg. \emph{totalResults}), in this way all the preceding steps will be performed automatically. Loading and handling alignments data require a lot of memory: if you have small datasets (up until 2.5Gb, ~58.5 million reads with a length of 50 bp) 4 Gb of RAM will be more than sufficient (peak memory usage, gotten via the ps command with the keyword rss, when comparing two samples of 2.5Gb: ~ 1.3Gb) and you can use the \emph{roarWrapper} script. For bigger datasets (eg. 20Gb, ~413 million reads with a length of 101bp) you will have to split your alignments in small chunks: \emph{roarWrapper\_chrBychr} works on single chromosomes. Please note than when working in this way FPKM values returned by \emph{fpkmResults} will be different from those obtained with the analysis performed using a single \Rclass{RoarDataset} on the whole alignments - in this case you can extract counts with \emph{countResults} and knit together real FPKM values afterwards, if you need them (\emph{roarWrapper\_chrBychr} does this). In the following sections we will show an example of this package usage on the \Rpackage{RNAseqData.HNRNPC.bam.chr14} package data and the annotation given in the package itself. \subsection{Creating a RoarDataset object} The analysis begins by creating a \Rclass{RoarDataset} object that holds the alignment data (4 HNRNPC ko samples and 3 control samples) and the annotation regarding 3'UTRs coordinates. <>= library(roar) library(rtracklayer) library(RNAseqData.HNRNPC.bam.chr14) gtf <- system.file("examples", "apa.gtf", package="roar") # HNRNPC ko data bamTreatment <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[seq(5,8)] # control (HeLa wt) bamControl <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[seq(1,3)] rds <- RoarDatasetFromFiles(bamTreatment, bamControl, gtf) @ \subsection{Obtaining counts} This is the first step in the Roar analysis: it counts reads overlapping with the PRE/POST portions defined in the given gtf/GRanges annotation. Reads of the given bam are accounted for with the following rules: \begin{enumerate} \item reads that align on only one of the given features are assigned to that feature, even if the overlap is not complete \item reads that align on both a PRE and a POST feature of the same gene (spanning reads) are assigned to the POST one, considering that they have clearly been obtained from the longer isoform \end{enumerate} If the alignments were derived from a strand-aware protocol it is possible to consider strandness when counting reads using the stranded=TRUE argument; in this case we use FALSE as long as we don't want to consider strandness. <>= rds <- countPrePost(rds, FALSE) @ \subsection{Computing roar} This is the second step in the Roar analysis: it computes the ratio of prevalence of the short and long isoforms for every gene in the treatment and control condition (m/M) and their ratio, roar, that indicates if there is a relative shortening-lengthening in a condition versus the other one (the choice of calling them ``treatment'' and ``control'' is simply a convention and reflects the fact that to calculate the roar the m/M of the treatment condition is used as the numerator). For details about the m/M calculations see section~\ref{appendix}. A roar > 1 for a given gene means that in the treatment condition that gene has an higher ratio of short vs long isoforms with respect to the control condition (and the opposite for roar < 1). Negative or NA m/M or roar occurs in not definite situations, such as counts equal to zero for PRE or POST portions. If for one of the conditions there is more than one sample, like in this example that has four and three samples for the two conditions, then calculations are performed on average counts. <>= rds <- computeRoars(rds) @ \subsection{Computing p-values} This is the third step in the Roar analyses: it applies a Fisher test comparing counts falling on the PRE and POST portion in the treatment and control conditions for every gene. If there are multiple samples for a condition every combination of comparisons between the samples lists is considered - for example in this case we have four and three samples, therefore it will calculate 12 p-values. <>= rds <- computePvals(rds) @ Computing all samples pairing is not the preferrable choice when a paired experimental design exists: in this case only the correct pairs between control and treatment samples should be compared with the Fisher test; then their p-values can be combined following the Fisher method (\cite{Fisher}) because we have different independent tests on the same null hypothesis. For these situations there is an alternative way to obtain the p-values: \emph{computePairedPvals}, which require the user to specify as arguments not only the \Rclass{RoarDataset} object but also two vectors containining the ordered samples for treatment and control - see the related man page for details. \subsection{Obtaining results} There are various functions aimed at extracting results from a \Rclass{RoarDataset}; the first one is \emph{totalResults}: <>= results <- totalResults(rds) @ It returns a dataframe with gene\_id as rownames and several columns: ``mM\_treatment'', ``mM\_control'', ``roar'', ``pval'' and a number of columns called ``pvalue\_X\_Y'' (when there are multiple samples for at least one condition, X refers to the treatment samples and Y to the control ones). The first two columns have the value of the ratio showing the relative abundance of the short isoform with respect to the long one in the treatment (or control) condition, the third one represents the roar and it's bigger than one when the shorter isoform is relatively more expressed in the treatment condition (and the other way around when it's smaller than one). The pval column stores the results of the Fisher tests when both conditions have a single sample, while the multiplication of all the possible combinations p-values for un-paired multiple samples analysis or the combined p-values following the Fisher method if \emph{computePairedPvals} has been used. In this case for example there will be columns ``pvalue\_1\_1'' up to ``pvalue\_4\_3'': the first number represents the number of the considered sample for the treatment condition, while the second for the control one (the samples are considered in the same order given in the lists passed to the \Rclass{RoarDataset} constructors). All other functions are simply helper functions: \emph{fpkmResults} adds columns "treatmentValue" and "controlValue" representing the level of expression of the relative gene in the two conditions, while \emph{countResults} puts in that columns the number of reads that were counted on the PRE portion of genes. \emph{standardFilter} get a parameter with a cutoff for the FPKM value and selects only genes with an expression higher than that and without any NA/Inf value for m/M or roar (deriving from not definite situations with counts equal to zero). \emph{pvalueFilter} further selects only genes with a Fisher test pvalue smaller than the given cutoff in the single sample case, while if more than one sample has been given for one condition it adds a column named "nUnderCutoff" that stores the number of p-values that are smaller than the given cutoff. It is worth noting that all the p-values are nominal with this method, while \emph{pvalueCorrectFilter} applies an user requested multiple correction procedure (after FPKM filtering) for single sample or paired (that is when \emph{computePairedPvals} has been used) analyses and filters with the given cutoff on these corrected value (for un-paired multiple samples design this does exactly the same as \emph{pvalueFilter}). For example in our case we can see how many genes have a pvalue under 0.05 in every comparisons in this way: <>= results_filtered <- pvalueFilter(rds, fpkmCutoff=-Inf, pvalCutoff=0.05) nrow(results_filtered[results_filtered$nUnderCutoff==12,]) @ Note that for FPKM we used -Inf as a cutoff because in this case we didn't want to filter out genes considering their expression values (and because we used only partial alignments deriving from a single chromosome, thus the usual FPKM cutoff are not directly applicable). Another thing to point out is that the FPKM values derive from the total reads falling over the genes given in the annotation, therefore when doing the analysis ``stepwise'' (eg. with roarWrapper\_chrByChr) their values will be different than those obtained performing the analysis all together. \emph{countResults} could be used in this situations to save single steps counts and then calculate ``whole'' FPKM values. \subsection{MultipleAPA analyses} In order to allow users to efficiently study polyadenylation for gene harbouring more than one alternative polyadenylation site with a single analysis we have developed the \Rclass{RoarDatasetMultipleAPA} class. A MultipleAPA analysis computes several roar values and p-values for each gene: one for every possible combinations of APA-canonical end of a gene (i.e. the end of its last exon). This is more efficient than performing several different ``standard'' roar analyses choosing the PRE and POST portions corresponding to different APAs because reads overlaps are computed only once. The gtf or \Rclass{Granges} used for this kind af analyses are different from the ones described previously: they should containing a suitable annotation of alternative APA sites and gene exonic structure. Minimal requirements are: attributes (metadata columns for \Rclass{Granges}) called "gene", "apa" and "type." APA should be single bases falling over one of the given genes and need to have the metadata column "type" equal to "apa" and the "apa" column with unambiguous ids and the corresponding gene id pasted together with an underscore. The "gene" metadata columns for these entries should not be initialized. All the studied gene exons need to be reported, in this case the metadata column "gene" should contain the gene id (the same one reported for each gene APAs) while "type" should be set to "gene" and "apa" to NA. All apa entries assigned to a gene should have coordinates that falls inside it and every gene that appears should contain at least one APA. Annotation files that follow these rules can be found at \url{http://94.91.177.36:6638/roar/}. From the user point of view the steps needed to perform a multipleAPA analysis are exactly the same of the standard one, the only difference lies in the output of methods that report results as long as in this case we have more than a single roar value for each gene. The details are reported in the corresponding man pages. \subsection{Genes with overlapping 3'UTRs} When two genes overlap reads that align on the shared coordinates could not be assigned to either of them in a straightforward way - in the single APA analyses we discard those reads while for technical reasons in the multiple APA analyses we perform the alignment separately for genes on the plus and minus strand, thus the reads that overlap genes on different strands will be counted for all of them, while those on genes on the same strand will be correctly discarded. In order to overcome this problem in the multipleAPA analyses it is possible to filter the gtf or \Rclass{Granges} annotation used before creating the \Rclass{RoarDatasetMultipleAPA} object removing all genes with overlapping portions. Note that for the single APA analyses roar alignment procedure could detect overlaps only between the given PRE and POST portions and therefore possible overlaps with other portions of genes are not considered - again it is possible to filter the used annotations beforehand. \section{Appendix}\label{appendix} \subsection{m/M calculations} To correctly evaluate the quantities ratios of the short and long isoforms counts of the reads falling over the two transcript portions (as previously defined: PRE is the portion common between the two isoforms, while POST is the one present only in the longer transcript) should be accounted for considering the fact that those falling over PRE could have been obtained from both isoforms while the others falling on POST derive exclusively from the longer one; another more trivial question that has to be addressed is that reads fall with larger frequencies on longer stretches of RNA. We can say that the total number of reads falling over a transcript in its entirety ($N$) derives from the relative abundance of the two isoforms and their potential of generating reads; that is: $ N = \epsilon_{M} M + \epsilon_{m} m$ where $m$ is the quantity of the short isoform, $M$ of the long one and $\epsilon$ identifies their efficiency in generating reads. Assuming the equiprobability of reads distribution (that is each nucleotide has the same probability of finding itself in a read) the efficiency in generating reads of the two isoforms is related to their lengths (with the addition of a constant of proportionality $k$):\vspace{0.4cm} \\ \begin{math} \epsilon_{M} = k (l_{PRE}+l_{POST}) \end{math} \vspace{0.4cm}\\ \begin{math} \epsilon_{m} = k (l_{PRE}) \end{math} \vspace{0.4cm}\\ Defining $l_{PRE}$ as the length of the PRE portion and $l_{POST}$ as the length of the POST we can now obtain the number of reads falling on the two portions as:\vspace{0.4cm}\\ \begin{math} \#r_{PRE} = \epsilon_{M} M (\frac{l_{PRE}}{l_{PRE}+l_{POST}}) + \epsilon_{m} m \end{math} \vspace{0.2cm}\\ and: \vspace{0.2cm}\\ \begin{math} \#r_{POST} = \epsilon_{M} M (\frac{l_{POST}}{l_{PRE}+l_{POST}}) \end{math} \vspace{0.5cm} These two equations reflect the fact that all the reads deriving from the short isoform ($\epsilon_{m} m$) fall on the PRE portion while the ones deriving from the long one are distributed over the PRE and POST portions depending on their lengths. We can now setup a system of equations aimed at obtaining the $m/M$ value in terms of the numbers of reads falling over the two portions and their lengths.\vspace{0.5cm} We start from: \begin{math} \left\{ \begin{array}{l} \#r_{PRE} = \epsilon_{M} M (\frac{l_{PRE}}{l_{PRE}+l_{POST}}) + \epsilon_{m} m \\ \#r_{POST} = \epsilon_{M} M (\frac{l_{POST}}{l_{PRE}+l_{POST}}) \end{array} \right. \end{math}\vspace{0.5cm}\\ Then with simple algebraic steps the system can be solved yielding the formula to obtain $m/M$ using only read counts and lengths:\vspace{0.5cm}\\ \begin{math} m/M = \frac{l_{POST} \#r_{PRE}}{l_{PRE} \#r_{POST}} -1 \end{math} \vspace{0.5cm} \\ As a simple emblematic case suppose that the PRE and POST portions have the same lengths and that the short and long isoforms are in perfect equilibrium (i.e. they are present in a cell in equal numbers). In this situation we will find on the PRE portion two times the number of reads falling on the POST one because half of them will derive from the long isoforms and the other half from the short ones. In this case the equation correcly gives us an $m/M$ equal to 1. In the previous discussion we have ignored reads falling across the PRE/POST boundaries. As long as they can derive only from the long isoform it is reasonable to assign them to $\#r_{POST}$. Portion lengths should be corrected to keep in consideration read lengths and the just cited assignment to POST of the spanning reads, therefore:\vspace{0.5cm}\\ \begin{math} l^{\prime}_{PRE} = l_{PRE} \end{math} \vspace{0.2cm} \\ \begin{math} l^{\prime}_{POST} = l_{POST} + readLength - 1 \end{math} \vspace{0.5cm} \\ % find a suitable relative measure XXX Normally we should have added readLength to both the lengths but in this case we do not expect reads to fall after the POST portion (that is the end of transcripts) and thus we only have to correct for the spanning reads. We do not subtract the same value from $l_{PRE}$ as long as in theory we could expect reads to fall at its 5' (i.e. reads falling across that exon and the previous one or those from still unspliced transcripts).\\ These corrected lengths are those used for the $m/M$ calculations. \bibliography{roar} % bib in the same dir as .Rnw \bibliographystyle{plain} \end{document}