RNA-seq: impact of RNA degradation on transcript quantification
© Gallego Romero et al.; licensee BioMed Central Ltd. 2014
Received: 23 May 2014
Accepted: 27 May 2014
Published: 30 May 2014
Skip to main content
© Gallego Romero et al.; licensee BioMed Central Ltd. 2014
Received: 23 May 2014
Accepted: 27 May 2014
Published: 30 May 2014
The use of low quality RNA samples in whole-genome gene expression profiling remains controversial. It is unclear if transcript degradation in low quality RNA samples occurs uniformly, in which case the effects of degradation can be corrected via data normalization, or whether different transcripts are degraded at different rates, potentially biasing measurements of expression levels. This concern has rendered the use of low quality RNA samples in whole-genome expression profiling problematic. Yet, low quality samples (for example, samples collected in the course of fieldwork) are at times the sole means of addressing specific questions.
We sought to quantify the impact of variation in RNA quality on estimates of gene expression levels based on RNA-seq data. To do so, we collected expression data from tissue samples that were allowed to decay for varying amounts of time prior to RNA extraction. The RNA samples we collected spanned the entire range of RNA Integrity Number (RIN) values (a metric commonly used to assess RNA quality). We observed widespread effects of RNA quality on measurements of gene expression levels, as well as a slight but significant loss of library complexity in more degraded samples.
While standard normalizations failed to account for the effects of degradation, we found that by explicitly controlling for the effects of RIN using a linear model framework we can correct for the majority of these effects. We conclude that in instances in which RIN and the effect of interest are not associated, this approach can help recover biologically meaningful signals in data from degraded RNA samples.
Degradation of RNA transcripts by the cellular machinery is a complex and highly regulated process. In live cells and tissues, the abundance of mRNA is tightly regulated, and transcripts are degraded at different rates by various mechanisms , partially in relation to their biological function [2–5]. In contrast, the fates of RNA transcripts in dying tissue, and the decay of isolated RNA are not part of normal cellular physiology and, therefore, are less likely to be tightly regulated. It remains largely unclear whether most transcript types decay at similar rates under such conditions or whether rates of RNA decay in dying tissues are associated with transcript-specific properties.
These questions are of great importance for studies that rely on sample collection in the field or in clinical settings (both from human populations as well as from other species), in which tissue samples often cannot immediately be stored in conditions that prevent RNA degradation. In these settings, extracted RNA is often partly degraded and may not faithfully represent in vivo gene expression levels. Sample storage in stabilizers like RNALater lessens this problem  but is not always feasible. Differences in RNA quality and sample handling could, therefore, confound subsequent analyses, especially if samples subjected to different amounts of degradation are naïvely compared against each other. The degree to which this confounder affects estimates of gene expression levels is not well understood.
There is also no consensus on the level of RNA decay that renders a sample unusable or on approaches to control for the effect of ex vivo processes in the analysis of gene expression data. Thus, while standardized RNA quality metrics such as the Degradometer  or the RNA Integrity Number (RIN; ), provide well-defined empirical methods to assess and compare sample quality, there is no widely accepted criterion for sample inclusion. For example, proposed thresholds for sample inclusion have varied between RIN values as high as 8  and as low as 3.95 . The recent Genotype-Tissue Expression (GTEx) project , for instance, reports both the number of total RNA samples they collected as well as the number of RNA samples with RIN scores higher than 6, presumably as a measure of the number of high quality samples in the study.
Broadly speaking, three approaches can be adopted to deal with RNA samples of variable quality. First, RNA samples with evidence of substantial degradation can be excluded from further study; this approach relies on establishing a cut-off value for ‘high quality’ versus ‘low quality’ samples and suffers from the current lack of consensus on what this cut-off should be. It also could exclude the possibility of utilizing unique and difficult to collect samples from remote locations or historical collections. Second, if investigators are willing to assume that all transcript types decay at a similar rate, variation in gene expression estimates due to differences in RNA integrity could be accounted for by applying standard normalization procedures. Third, if different transcripts decay at different rates, and if these rates are consistent across samples for a given level of RNA degradation – for example, a given RIN value – a model that explicitly incorporates measured, sample-specific, degradation levels could be applied to gene expression data to correct for the confounding effects of degradation.
To date, most studies apply a combination of the first two approaches: an application of an arbitrary RNA quality cutoff (typically based on RIN score), followed by standard normalization of the data, which assumes that RNA samples at any RIN value higher than the chosen cutoff are not subjected to transcript-specific decay rates. However, current work on the effects of RNA decay has not yet provided clear guidelines with respect to these approaches. In addition, nearly all published work that focuses on RNA stability in tissues following cell death and/or sample isolation predates, or does not employ, high throughput sequencing technologies. These studies broadly suggest that both the quantity and quality of recovered RNA from tissues can be affected by acute pre-mortem stressors, such as pyrexia or prolonged hypoxia [12–14], and by the time to sample preservation and RNA extraction. The quantity and quality of recovered RNA are strongly dependent on the type of tissue studied , even when sampling from the same individual [16, 17]. These differences in yield across tissues have resulted in a wide range of recommendations for an acceptable post-mortem interval for extracting usable, high-quality RNA, ranging from as little as 10 minutes  to upwards of 48 hours , depending on tissue source and preservation conditions.
Similarly, studies examining changes in the relative abundance of specific transcripts as a result of ex vivo RNA decay have reached somewhat contradictory recommendations. Some of this conflict may be attributable to methodological differences. Studies that focused on small numbers of genes assayed through quantitative PCR consistently report little to no effect of variation in RNA quality on gene expression estimates [6, 19–22]. Conversely, microarray-based studies have repeatedly reported significant effects of variation of RNA quality on gene expression estimates, even after applying standard normalization approaches. Increasing the time from tissue harvesting to RNA extraction or cryopreservation from 0 to only 40 or 60 minutes, for example, significantly affected expression profiles in roughly 70% of surveyed genes in an experiment on human colon cancer tissues . Likewise, a substantial fraction of genes in peripheral blood mononuclear cells (PBMCs) appears to be sensitive to ex vivo incubation . Other microarray-based studies have reached similar conclusions, both in samples from humans [15, 16, 22, 23] and other organisms , and have urged caution when analyzing RNA samples with medium or low RIN scores, although the definition of an acceptable RNA quality threshold remains elusive.
To examine the effects of RNA degradation in a setting relevant to field study sample collection, we sequenced RNA extracted from PBMC samples that were stored unprocessed at room temperature for different time periods, up to 84 hours. We collected RNA decay time-course data spanning almost the entire RIN quality scale and examined relative gene-specific degradation rates through RNA sequencing. Due to the high sensitivity and resolution of high-throughput RNA sequencing, our data provide an unprecedentedly detailed picture of the dynamics of RNA degradation in stressed, ex vivo cells. Based on our results, we develop specific recommendations for accounting for these effects in gene expression studies.
We extracted RNA from 32 aliquots of PBMC samples from four individuals. The PBMC samples were stored at room temperature for 0 hours, 12 hours, 24 hours, 36 hours, 48 hours, 60 hours, 72 hours and 84 hours prior to RNA extraction. As expected, time to extraction significantly affected the RNA quality (P <10−11), with mean RIN = 9.3 at 0 hours and 3.8 at 84 hours [see Additional file 1: Table S1]. Based on the RIN values we chose to focus on 20 samples from five time points (0 hours, 12 hours, 24 hours, 48 hours and 84 hours) that spanned the entire scale of RNA quality. We generated poly-A-enriched RNA sequencing libraries from the 20 samples using a standard RNA sequencing library preparation protocol (see ). We added a spike-in of non-human control RNA to each sample, which allowed us to confirm the effects of RNA degradation on the RNA sequencing results (see Methods for more details). Following sequencing, we randomly subsampled all libraries to a depth of 12,129,475 reads, the lowest number of reads/library observed in the data. We used BWA 0.6.3 to map reads, calculated reads per kilobase transcript per million (RPKM), and normalized the data using a standard quantile normalization approach (for example, as in ). We observed that sample RIN is associated with both the number of uniquely mapped reads (analysis of variance (ANOVA) P <10−3) and the number of reads mapped to genes (P <10−3; Additional file 2: Figure S1), with high RIN samples having greater numbers of both. Furthermore, the proportion of exogenous spike-in reads increases significantly as RIN decreases (P <10−10), as expected given degradation-driven loss of intact human transcripts in poor quality samples. Sequence reads from individual 2 were poorly mapped, especially in the later time-points (see Methods and Additional file 2: Figure S1); we thus excluded the data from all samples from this individual in subsequent analysis.
Genes observed in all individuals until or after a particular time point
Mean RPKM when seen
Mean RPKM when seen
Although we might expect RNA degradation in decaying cells to be a random process, gene ontology (GO) analysis identified 118 and 293 significantly overrepresented categories among slowly and rapidly degraded genes, respectively (FDR = 5%; Additional file 12: Tables S4 and Additional file 13: Table S5). We present the functional enrichment results only as an indication that the rate of transcript decay is not random. These observations are robust to different normalization approaches, to the inclusion of RIN as a covariate in our linear model, and to fitting slopes using RIN instead of time-points. Limiting our analyses to the 1,000 bp closest to the 3′ end of transcripts also yields similar results.
We also sought to investigate whether targets of slow, fast, or average degradation differ meaningfully in terms of broad biological categories. As expected given our poly-A enrichment strategy, most transcripts in our data originate from intact protein-coding genes, but we also observed four other biotypes represented by more than 100 distinct transcripts. The distribution of these biotypes across rapidly and slowly degraded transcripts is not random, with a significant enrichment of pseudogenes among transcripts that degrade slowly (P = 0.015), and an enrichment of intact protein-coding genes among the rapidly degraded transcripts (P <10−16, Figure 4E).
Ultimately, the goal of most RNA sequencing studies is to estimate variation in gene expression levels or to identify genes that are differentially expressed between conditions, individuals or states. We thus considered the effects of RNA quality on measures of relative gene expression levels between time-points and on overall estimates of inter-individual variation in gene expression.
Number of identified DE genes
GLM: reads approximate time point
GLM + RIN
Regress RIN, GLM
0 h versus 12 h
0 h versus 24 h
0 h versus 48 h
0 h versus 84 h
GLM: reads approximate individual b
GLM + RIN
Regress RIN, GLM
Ind 1 versus Ind 3
Ind 1 versus Ind 4
Ind 3 versus Ind 4
As expected, when we include RIN as a covariate in the model the number of differentially expressed genes across time-points is drastically reduced (fewer than 50 genes are classified as differentially expressed between 0 hours and any other time-point; Table 2). These observations confirm that RIN is a robust indicator of degradation levels. Without accounting for RIN, the effect of variation in RNA quality on our data is overwhelming. To understand these effects better, we explored whether accounting for variation in RIN increased the power to detect other sources of (biologically relevant) variation in RNA-seq data, such as the variation in gene expression between individuals. We also investigated several alternative approaches for controlling for variation in RNA quality.
To recover this signal, we tested two approaches for explicitly accounting for RIN when estimating differential gene expression across individuals: (1) incorporating RIN as a covariate in our GLM; and (2) analyzing the residuals of gene expression levels after first regressing out RIN from the normalized gene expression data (Table 2). Both approaches result in the identification of many more genes as differentially expressed between individuals (401 to 573 when incorporating RIN directly into our GLM, 190 to 299 when testing for differential expression using residuals; Table 2). We also repeated the pairwise correlation analysis using the same 1,410 most variable genes identified above, but this time we used the residuals after regressing the effect of RIN from the data. The residuals cluster well by individual throughout the entire time course experiment, regardless of RNA quality (Figure 5B).
DE genes across pairs of individuals and overlap with top 10% most variable genes at 0 hours
GLM, individual + RIN
Regress RIN, GLM individual
Number DE genes
Number DE genes
Number DE genes
Ind 1 vs Ind 3
Ind 1 vs Ind 4
Ind 3 vs Ind 4
Our observations indicate that the effects of RNA degradation following death or tissue isolation are pervasive and can rapidly obscure inter-individual differences in gene expression. Yet, we also found that by using RNAseq nearly all genes observed at our first time-point could still be detected even in severely degraded RNA samples, although the estimated relative expression levels were drastically affected by degradation. Although postmortem RNA degradation is considered a non-regulated process, some of the traditional predictors of regulated RNA decay rates in the cell are also associated with variation in RNA quality in our data. For example, longer protein coding regions and 3′ UTRs are correlated with more rapid degradation, similar to previously reported trends [5, 28, 29]. Total transcript length, however, which is a significant predictor of regulated RNA decay in the cell, is not associated with variation in degradation rates in our data.
We confirmed previous observations of decreasing data quality as time from tissue extraction to RNA isolation increased [see Additional file 2: Figure S1], both with respect to the number of high quality reads we were able to generate from our sequencing libraries and library complexity. While increased time to RNA extraction did not generally result in the complete loss of transcripts (less than 8% of transcripts are lost), the relative expression levels of many transcripts were drastically altered over the time-course experiment, with 61% of genes classified as differentially expressed between 0 hours (mean RIN of 9.3) and 84 hours (mean RIN of 3.78). This proportion of differentially expressed genes is in line with previous reports of the effects of warm ischemia on human gene expression in tumor biopsies, as assessed using microarrays [20, 22]. The potential of RNA degradation to skew measurements of gene expression levels and obscure biologically meaningful signals is, therefore, apparent. If there are systematic differences in RNA quality between two classes of samples being compared, we predict that the effect of RNA quality on relative estimates of gene expression levels would be responsible for much of the signal in the data. Furthermore, as degradation rate is to some degree associated with biological function [see Additional file 11: Tables S3 and Additional file 12: Table S4], it has the potential to confound naïve comparisons of functional annotations as well.
However, the marked effects of RNA degradation on the relative expression level of most genes can be corrected, to a large degree, using relatively simple statistical methods. Indeed, we found that the inclusion of RIN in our model was sufficient to account for much of the effect of degradation and allowed us to identify a reasonable number of differentially expressed genes between pairs of individuals in our data. These were not spurious signals generated by our approach; they recapitulated observations made at 0 hours (when RNA quality was excellent), but were originally dwarfed by the magnitude of degradation-driven expression changes in the uncorrected data. A similar approach – taking into account variation in RIN – has been previously proposed for the analysis of RTq-PCR data abundance . Nevertheless, our observations suggest that some of the effects of transcriptional degradation in ex vivo samples cannot be corrected. Library complexity decreases somewhat with lower RNA quality, and some genes (approximately 5%) can no longer be detected at the later time-points. Based on our data we conclude that these effects cannot be corrected by simply sequencing more degraded libraries to a greater depth.
In a study similar to our own, Opitz et al.  subjected extracted RNA samples from three advanced human rectal cancer biopsies to degradation through increasingly longer incubation at 60˚C and then considered the evidence of time-point/RIN–driven degradation using microarray data. The RIN values spanned by their data mirror values in ours, but the results do not. In contrast to the large RIN-associated effects we observed, Opitz et al. reported that of 41,000 tested probe-level 2data points only 275 demonstrated significant degradation effects, with inter-individual differences being the predominant signal in the data. Assuming that differences in the platforms used (microarrays and RNAseq) are not the reason for this discrepancy, one possible explanation for this stark difference between the studies is that lower RIN scores as a result of degradation of extracted RNA samples (Opitz et al.) may reflect substantially different properties than lower RIN scores that are the result of degradation of RNA in decaying cells (our study). Based on the observations of Opitz et al. we hypothesize that degradation rates of isolated RNA may be mostly linear and uniform; thus, the degradation effects can be accounted for by employing standard normalization approaches. In contrast, degradation rates of RNA in a dying tissue sample, a situation that mirrors more closely conditions likely to be faced by investigators in clinical or field settings, is not uniform across transcripts. Because these differences cannot be neglected in downstream analyses, knowledge of the context in which degradation occurs is, therefore, crucial.
Our observations suggest that actively mediated degradation of transcripts may occur during necrosis; namely, degradation of RNA in a dying tissue may not be a completely random process. Biologically mediated degradation, whether actively driven by the cell’s decay machinery , or simply the consequence of the leakage of RNases into cells as membranes are disrupted, is different from the heat-driven degradation of naked RNA, which in turn is likely to be different from the degradation caused by continued freeze-thaw cycles . It is likely that in a dying tissue, most degradation is initially biologically mediated and directed towards specific classes of transcripts, but as the cellular environment continues to deteriorate, the relative importance of stochastic degradation may increase such that at later time-points degradation becomes increasingly uncoupled from biological function.
Additionally, the increased resolution of RNA sequencing relative to other platforms used to assay gene expression levels  is both a hindrance and a boon in this situation, allowing for detection of subtler differences than ever before, but also warranting greater caution when analyzing samples of differing quality.
Previous studies [8–10, 23, 32–34] have sought to provide an RNA degradation threshold below which in-depth analysis of RNA is not recommended. However, these studies have reached conflicting conclusions. Our data suggest that if a simple cut-off value is to be used, a conservative cut-off in the context of RNA degradation in dying tissue samples lies between 7.9 and 6.4, the mean RIN scores associated with 12 hours and 24 hours in our time course experiment, respectively. We observed few differences in measurements of gene expression between 0 hours and 12 hours, as evidenced by the low number of genes identified as differentially expressed between the two time-points. Thus, it may be tempting to conclude that so long as all samples in any particular study have roughly similar RINs explicit correction is not necessary. However, when we test for differential expression between other close time-points we identify 3,020 genes as differentially expressed between 48 hours and 84 hours (difference in mean RIN = 1.3), and 5,293 between 24 hours and 48 hours (difference in mean RIN = 1). It is clear that measurements of gene expression are extremely sensitive to starting sample quality.
Our observations indicate that useful data can be collected using RNA sequencing even from highly degraded samples. As long as RIN scores are not associated with the effect of interest in the study (namely, different classes of samples in the study are not associated with different distributions of RIN scores), accounting for RIN scores explicitly can be an effective approach. In our study, we were able to identify differently expressed genes between individuals even when RNA samples with RIN scores around 4 were included. Excluding the samples with RIN values lower than 6.4 in our study would have resulted in a less powerful design than including these samples and globally correcting for RIN values. Given these results, we believe that under most circumstances, the most effective approach may be to include all samples regardless of quality, and explicitly model a measure of RNA quality in the analysis.
We obtained Buffy coat samples from four adult Caucasian males from Research Blood Components LLC (Boston, MA, USA) and separated PBMCs through a standard Ficoll gradient purification. Each sample was split into aliquots of four million live cells and resuspended in 200 uL of PBS. Cells were kept at room temperature and aliquots from each sample lysed every twelve hours by addition of 700 uL of RLT buffer (Qiagen, Valencia, CA, USA) with beta-mercaptoethanol (Sigma-Aldrich, St Louis, MO, USA) added at 10 uL BME/1 mL RLT according to the manufacturer’s instructions. Lysed cells were immediately frozen and not thawed until RNA extraction.
RNA was extracted using the Qiagen RNeasy kit. Extracted RNA quality was assessed with a BioAnalyzer (Agilent Technologies, Wilmington, DE, USA). From these results we selected five time-points – 0 hours, 12 hours, 24 hours, 48 hours and 84 hours – that encompassed a large stretch of the degradation spectrum. We then prepared poly-A-enriched RNA sequencing libraries for all 20 individual/time-point combinations according to a previously published protocol , using 1.5 μg of total RNA per library in all instances. In all instances, we added 15 ng (1%) of an exogenous RNA spike-in during library preparation, composed of equal parts Caenorhabditis elegans, Drosophila melanogaster and Danio rerio total RNA. Samples were multiplexed and sequenced on four lanes (two per library preparation strategy) of an Illumina HiSeq2000 using standard protocols and reagents. Reads were 50 bp in length. All generated reads have been deposited into the Sequence Read Archive (SRA) under accession numbers SAMN02769865-SAMN02769884.
Data were combined across lanes and data for all libraries were randomly subsampled to the lowest observed number of reads, 12,129,475. Reads were independently mapped to the human (hg19), D. rerio (danRer7), D. melanogaster (dm3), and C. elegans (ce10) genomes using BWA 0.6.2 . All reference genomes were obtained from the UCSC Genome Browser . Only reads that mapped exclusively to a single site in the human genome with one or zero mismatches were retained for downstream analyses. Following mapping, we removed all reads that mapped to more than one genome. At this point we also discarded one individual – individual number 2 - due to low mappability and read quality in the later time points [see Additional file 2: Figure S1]. We also mapped all reads using TopHat 2.0.8 and the same quality thresholds and filtering steps.
We calculated RPKM  for all ENSEMBL v71  human genes in our data. Genes with multiple transcripts were collapsed into a single transcript containing all exons of the gene; where multiple exons of different size overlapped the same genomic region, the entire region was kept. We discarded all exonic regions transcribed as part of more than one gene. Additionally, we quantile-normalized both RPKM and read count-level data across individuals using the lumiN function in the Bioconductor  package lumi , which controls for, and dampens, technical sampling variance in highly expressed genes. Read counts were log2 transformed prior to quantile normalization to generate a normal distribution; analyses were carried out on subsequently untransformed counts.
All statistical analyses were carried out using R 2.15.2.
where y(t) is the mRNA abundance of a given gene at time t (in quantile-normalized RPKM), B0 is the abundance at the initial time-point, and k the decay rate, with the variance term ϵ being normally distributed. Data from all three individuals were considered simultaneously; that is, we obtained a single decay constant for each gene across all three individuals. To control for the high FDR of expressed genes at low expression levels, all RPKM observations <0.3 were discarded for all subsequent analyses, as in .
Length and per-transcript %GC content were calculated using BEDTools (version 2.16.2 ), using the same gene models described above. Biotype as well as 5′ and 3′ UTR length were retrieved from ENSEMBL v71. In those instances where there are multiple UTRs associated with the same gene, we used the median UTR length for each gene in all calculations.
Differentially expressed genes were identified using the R package edgeR , utilizing a GLM framework with time, individual ID and sample RIN as covariates, as described above. Only those genes with a minimum observation of one mapped read across all individuals at all time-points were included. Instead of quantile normalization as described above, all data were normalized using trimmed mean of M values normalization (TMM, ), which corrects for the observed differences in informative reads between sequencing libraries. Inter-individual variance estimates were generated after variance stabilization of read counts using the predFC function in edgeR.
Downstream gene enrichment analyses were carried out using the R package topGO , using the ‘classic’ algorithm and a minimum node size of five. All significance values given in the text have been corrected to an FDR of 5% or 1%, using the qvalue method of . In all cases, the background data set included all 14,094 genes with complete observations.
We would like to thank Rachel A. Myers as well as members of the Gilad lab for helpful comments and discussions. This work was supported by NIH grant HG006123 to YG. IGR was supported by a Sir Henry Wellcome Postdoctoral Fellowship. AAP was supported by an American Heart Foundation Predoctoral Fellowship. JT was supported by a Chicago Fellows Postdoctoral Fellowship.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.