Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences

Introduction

Quantification and comparison of isoform- or gene-level expression based on high throughput sequencing reads from cDNA (RNA-seq) are arguably among the most common tasks in modern computational molecular biology. Currently, one of the most widely used approaches amounts to defining the genomic locations of a set of non-overlapping targets (typically, genes) and using the number of aligned reads overlapping a target as a measure of its abundance, or expression level. Several software packages have been developed for performing such “simple” counting (e.g., featureCounts¹ and HTSeq-count²). More recently, the field has seen a surge in methods aimed at quantifying the abundances of individual transcripts (e.g., Cufflinks³, RSEM⁴, BitSeq⁵, kallisto⁶ and Salmon⁷). These methods provide higher resolution than simple counting, and by circumventing the computationally costly read alignment step, some (notably, kallisto and Salmon) are also considerably faster. However, isoform quantification is more complex than the simple counting, due to the high degree of overlap among transcripts. Currently, there is no consensus regarding the optimal resolution or method for quantification and downstream analysis of transcriptomic output.

Another point of debate is the unit in which abundances are given. The traditional R/FPKM^8,9 (reads/fragments per kilobase per million reads) have been largely superseded by the TPM¹⁰ (transcripts per million), since the latter is more consistent across libraries. Regardless, all these units attempt to “correct for” sequencing depth and feature length and thus do not reflect the influence of these on quantification uncertainty. In order to account for these aspects, most statistical tools for analysis of RNA-seq data operate instead on the count scale. Most of these tools were designed to be applied to simple read counts, and the degree to which their performance is affected by using fractional estimated counts resulting from portioning reads aligning to multiple transcripts is still an open question. The fact that the most common sequencing protocols provide reads that are much shorter than the average transcript implies that the observed read counts depend on a transcript’s length as well as its abundance; thus, simple counts are arguably less accurate measures than TPMs of the true abundance of RNA molecules from given genes. The use of gene counts as input to statistical tools typically assumes that the length of the expressed part of a gene does not change across samples and thus its impact can be ignored for differential analysis.

In the analysis of transcriptomic data, as for any other application, it is of utmost importance that the question of interest is precisely defined before a computational approach is selected. Often, the interest lies in comparing the transcriptional output between different conditions, and most RNA-seq studies can be classified as either: 1) differential gene expression (DGE) studies, where the overall transcriptional output of each gene is compared between conditions; 2) differential transcript/exon usage (DTU/DEU) studies, where the composition of a gene’s isoform abundance spectrum is compared between conditions, or 3) differential transcript expression (DTE) studies, where the interest lies in whether individual transcripts show differential expression between conditions. DTE analysis results can be represented on the individual transcript level, or aggregated to the gene level, e.g., by evaluating whether at least one of the isoforms shows evidence of differential abundance.

In this report, we make and give evidence for three claims: 1) gene-level estimation is considerably more accurate than transcript-level; 2) regardless of the level at which abundance estimation is done, inferences at the gene level are appealing in terms of robustness, statistical performance and interpretation; 3) taking advantage of transcript-level abundance estimates when defining or analyzing gene-level abundances leads to improved DGE results compared to simple counting for genes exhibiting DTU. The magnitude of the effect in a given data set thus depends on the extent of DTU, and the global impact is relatively small in several real data sets analyzed in this study.

To facilitate a broad range of analysis choices, depending on the biological question of interest, we provide an R/Bioconductor package, tximport, to import transcript lengths and abundance estimates from several popular quantification packages and export (estimated) count matrices and, optionally, average transcript length correction terms (i.e., offsets) that can be used as inputs to common statistical engines, such as DESeq2¹¹, edgeR¹² and limma¹³.

Data and methods

Throughout this manuscript, we utilize two simulated data sets and four experimental data sets (Bottomly¹⁴ [Data set 3], GSE64570¹⁵ [Data set 4], GSE69244¹⁶ [Data set 5], GSE72165¹⁷ [Data set 6], see Supplementary File 1 for further details) for illustration. Details on the data generation and full records of the analyses are provided in the data sets and Supplementary File 1. The first simulated data set (sim1; Data set 1) is the synthetic human data set from Soneson et al.¹⁸, comprising 20,410 genes and 145,342 transcripts and is available from ArrayExpress (accession E-MTAB-3766). This data set consists of three biological replicates from each of two simulated conditions, and differential isoform usage was introduced for 1,000 genes by swapping the relative expression levels of the two most dominant isoforms between conditions. For each gene in this data set, the total transcriptional output is the same in the two conditions (i.e., no overall DGE); it is worth noting that this is an extreme situation, but provides a useful test set for contrasting DGE, DTU and DTE. The second simulated data set (sim2; Data set 2) is a synthetic data set comprising the 3,858 genes and 15,677 transcripts from the human chromosome 1. It is available from ArrayExpress with accession E-MTAB-4119. Also here, we simulated two conditions with three biological replicates each. For this data set, we simulated both overall DGE, where all transcripts of the affected gene showed the same fold change between the conditions (420 genes), differential transcript usage (DTU), where the total transcriptional output was kept constant but the relative contribution from the transcripts changed (420 genes) and differential transcript expression (DTE), where the expression of 10% of the transcripts of each affected gene was modified (422 genes, 528 transcripts). The three sets of modified genes were disjoint. Again, this synthetic data set represents an extreme situation compared to most real data sets, but provides a useful test case to identify underlying causes of differences between results from various analysis pipelines.

In addition to Data set 1–Data set 6, which contain all code for reproducing our analyses, further method descriptions are given in Supplementary File 1.

Gene abundance estimates are more accurate than transcript abundance estimates

To evaluate the accuracy of abundance estimation with transcript and gene resolution, we used the quasi-mapping mode of Salmon⁷ (v0.5.1) to estimate the TPM for each transcript in each of the data sets. Gene-level TPM estimates, representing the overall transcriptional output of each gene, were obtained by summing the corresponding transcript-level TPM estimates. For the two simulated data sets, the true underlying TPM of each feature was known and we could thus evaluate the accuracy of the estimates. Unsurprisingly, gene-level estimates were more accurate than transcript-level estimates (Figure 1A, Supplementary Figures 4,5). We also derived TPM estimates from simple gene-level counts obtained from traditional alignment of the reads to the genome using STAR followed by counting with featureCounts, by dividing the read count for each gene with a reasonable measure of the length of the gene (the length of the union of its exons) and the total number of mapped reads, and scaling the estimates to sum to 1 million. The simple count estimates showed a lower correlation with the true TPMs than the Salmon estimates, in line with previous observations¹⁹. It is worth noting that we are comparing entire (typical) workflows, and that differences may also occur if the set of reads that STAR is able to align to the genome is not identical to the set of reads that are contributing to the abundance estimation of Salmon. However, due to the large fraction of aligned reads and the high mapping rate with Salmon (both exceeding 99.8%, more than 95% of the reads were subsequently unambiguously assigned to genes by featureCounts), we do not expect this to have a major impact on the results shown in Figure 1A.

Figure 1 (sim2).

A: Accuracy of gene- and transcript-level TPM estimates from Salmon and scaled FPKM estimates derived from simple counts from featureCounts, in one of the simulated samples (sampleA1). Spearman correlations are indicated in the respective panels. Top row: using the complete annotation. Bottom row: using an incomplete annotation, with 20% of the transcripts randomly removed. Gene-level estimates are more accurate than transcript-level estimates. Gene-level estimates from Salmon are more accurate than those from featureCounts. B: Distribution of the coefficients of variation of gene- and transcript-level abundance estimates from Salmon, calculated across 30 bootstrap samples of one of the simulated samples (sampleA1). Gene-level estimates are less variable than transcript-level estimates. C: An example of unidentifiable transcript-level estimates, as uneven coverage does not cover the critical regions that would determine the amount that each transcript is expressed, while gene-level estimation is still possible.

Gene-level estimates derived from both simple counts and Salmon tended to show a high degree of robustness against incompleteness of the annotation catalog, as evidenced from estimation errors after first removing (at random) 20% of the transcripts (Figure 1A, see also Supplementary File 1); in contrast, Salmon’s transcript estimate accuracies deteriorated. To further compare the merits of genome alignment-based vs alignment-free quantification, especially in their handling of multi-mapping reads, we investigated the accuracy of the abundance estimates within sets of paralogous genes (Supplementary Figures 1–3). Also here, Salmon provided more consistently accurate estimates than STAR+featureCounts. From the bootstrap estimates generated by Salmon, we also estimated the coefficient of variation of the abundance estimates. The gene-level estimates showed considerably lower variability than the transcript-level estimates in both simulated and experimental data (Figure 1B, Supplementary Figures 6,7). Taken together, these observations suggest that the gene-level estimates are more accurate than transcript-level estimates and therefore potentially allow a more accurate and stable statistical analysis. A further argument in favor of gene-level analysis is the unidentifiability of transcript expression that can result from uneven coverage caused by underlying technical biases. While some extent of coverage variability might be alleviated by corrections for sequence- or position-specific biases²⁰, there remain cases where transcript expression cannot be inferred from data (Figure 1C). Intermediate approaches, grouping together “indistinguishable” features are also conceiveable²¹, but not yet standard practice.

DTE is more powerful and easier to interpret on gene level than for individual transcripts

DTE is concerned with inference of changes in abundance at transcript resolution, and thus invokes a statistical test for each transcript. We argue that this can lead to several complications: the first is conceptual, since the rows (transcripts) in the result table will in many cases not be interpreted independently, since the researcher is often interested in comparing the results for transcripts from the same gene locus, and the second one is more technical, since the number of transcripts is considerably larger than the number of genes, which could lead to lower power due to the portioning of the total set of reads across a larger number of features and a potentially higher multiple testing penalty. We tested for DTE on the simulated data by applying edgeR¹² to the transcript counts obtained from Salmon (the application of count models to estimated counts is discussed in the next Section), and represented the results as transcript-level p-values or aggregated these to the gene level by using the perGeneQValue function from the DEXSeq²² R package. Note that the transcript-level DTE test assesses the null hypothesis that individual transcripts do not change their expression, whereas the gene-level DTE test assesses the null hypothesis that all transcripts from a given gene exhibit no change in expression. Framing the DTE question at the gene level results in higher power, without sacrificing false discovery rate control (Figure 2A). This is not surprising given the different null hypotheses, and, in fact, for many of the genes detected as true positives with the gene-level test, only a subset of the truly changing transcripts were detected (Supplementary Figure 8). We note that this type of gene-level aggregation may favor genes in which one transcript shows strong changes, and that other approaches to increase power against specific alternatives are conceivable, e.g., capitalizing on the rich collection of methods for gene set analysis.

Figure 2 (sim2).

A: DTE detection performance on transcript- and gene-level, using edgeR applied to transcript-level estimated counts from Salmon. The statistical analysis was performed on transcript level and aggregated for each gene using the perGeneQValue function from the DEXSeq R package; aggregated results show higher detection power. The curves trace out the observed FDR and TPR for each significance cutoff value. The three circles mark the performance at adjusted p-value cutoffs of 0.01, 0.05 and 0.1. B: Schematic illustration of different ways in which differential transcript expression (DTE) can arise, in terms of absence or presence of differential gene expression (DGE) and differential transcript usage (DTU).

While DTE analysis is more suitable than DGE analysis for detecting genes with changes in absolute or relative isoform expression but no or only minor change in overall output (Supplementary Figure 9), we argue that even gene-level DTE results may suffer from lack of interpretability. DTE can manifest in several different ways, as an overall differential expression of the gene or differential relative usage of its transcripts, or a combination of the two (Figure 2B). We argue that the biological question of interest is in many cases more readily interpretable as a combination of DGE and DTU, rather than DTE. It has been our experience that results reported at the transcript level are still often cast to the gene level (i.e., given a differentially expressed transcript, researchers want to know whether also other isoforms of the gene are changing), suggesting that asking two specific gene-level questions (Is the overall abundance changing? Are the isoform abundances changing proportionally?) trumps the interpretability of one broad question addressing the transcript abundances (Are there changes in any of the isoform expression levels?), despite the increased need for multiple testing correction associated with performing two tests for each gene rather than one. There are of course also situations when a transcript-centric approach provides superior interpretability, for example in targeted experiments where specific isoforms are expected to change due to an administered treatment.

Incorporating transcript-level estimates leads to more accurate DGE results

DGE (i.e., testing for changes in the overall transcriptional output of a gene) is typically performed by applying a count-based inference method from statistical packages such as edgeR¹² or DESeq2¹¹ to gene counts obtained by read counting software such as featureCounts¹, HTSeq-count² or functions from the GenomicAlignments²³ R package. A lot has been written about how simple counting approaches are prone to give erroneous results for genes with changes in relative isoform usage, due to the direct dependence of the observed read count on the transcript length²⁴, and alternatives, such as Cuffdiff²⁴, which utilizes estimated transcript abundances, have been proposed. However, the extent of the problem in real data has not been thoroughly investigated. Here, we show that taking advantage of transcript-resolution estimates (e.g., obtained by Salmon) in count-based inference methods can lead to improved DGE results. We propose two alternative ways of integrating transcript abundance estimates into the DGE pipeline: to define an “artificial” count matrix, or to calculate offsets that can be used in the statistical modeling of the observed gene counts from, e.g., featureCounts. Both approaches are implemented in the accompanying tximport R package (available from http://bioconductor.org/packages/tximport).

For the DGE analyses, we defined three different gene-level count matrices for each data set (see also Supplementary File 1): 1) using featureCounts from the Rsubread¹ R package (denoted featureCounts below), 2) summing the estimated transcript counts from Salmon within genes (simplesum), 3) summing the estimated transcript TPMs from Salmon within genes, and multiplying with the total library size in millions (scaledTPM). We note that the scaledTPM values are artificial values, transforming underlying abundance measures to the “count scale” to incorporate the information provided by the sequencing depth. We further used the effective transcript lengths and estimated TPMs from Salmon to define average transcript lengths for each gene and each sample (normalization factors) as described in the Supplementary material, to be used as offsets for edgeR and DESeq2 when analyzing the featureCounts and simplesum count matrices (featureCounts_avetxl and simplesum_avetxl).

Overall, the counts obtained by all methods were highly correlated (Supplementary Figures 10–12), which is not surprising since any differences are likely to affect a relatively small subset of the genes. In general, the simplesum and featureCounts matrices led to similar conclusions in all considered data sets, even though there are differences between the two approaches in terms of how multi-mapping reads and reads partly overlapping intronic regions are handled²⁵. Previous studies have also shown that some loss of sensitivity for certain genes may be encountered from discarding multi-mapping fragments, which may be recovered through the use of transcript abundance estimators such as Salmon²¹. The concordance between simplesum and featureCounts results also suggests that statistical methods based on the Negative Binomial assumption are applicable also to summarized, gene-level estimated counts, which is further supported by the similarity between the p-value histograms as well as the mean-variance relationships observed with the three types of count matrices (Supplementary Figures 13–18).

Accounting for the potentially varying average transcript length across samples when performing DGE, either in the definition of the count matrix (scaledTPM) or by defining offsets (featureCounts_avetxl, simplesum_avetxl), led to considerably improved false discovery rate (FDR) control compared to using the observed featureCounts or aggregated Salmon counts directly (Figure 3A, Table 1). It is important to note that this improvement is entirely attributable to an improved handling of genes with changes in isoform composition between the conditions (Figure 3B, Supplementary Figure 19), that we purposely introduced strong signals in the simulated data set in order to pinpoint these underlying causes, and that the overall effect in a real data set will depend on the extent to which considerable DTU is present. Experiments on various real data sets (Supplementary Figure 20) show only small differences in the collections of significant genes found with the simplesum and simplesum_avetxl approaches, suggesting that the extent of the problem in many real data sets is limited, and that most findings obtained with simple counting are not induced by counting artifacts. Further support for this conclusion is shown in Figure 4 (see also Supplementary Figures 21–23 and Supplementary Table 1), where log-fold change estimates from edgeR, based on the simplesum and scaledTPM matrices, are contrasted. For the genes with induced DTU in the sim2 data set, log-fold changes based on the simplesum matrix are overestimated, as expected. However, this effect is almost absent in all the real data sets, again highlighting the extreme nature of our simulated data and suggesting that the effect of using different count matrices is considerably smaller for many real data sets. Table 1 further suggests that the lack of error control for simplesum and featureCounts matrices is more pronounced when there is a large difference in length between the differentially used isoforms. In the group with smallest length difference, where the longer differentially used isoform is less than 34% longer than the shorter one, all approaches controlled the type I error satisfactorily. It is worth noting that among all human transcript pairs in which both transcripts belong to the same gene, the median length ratio is 1.85, and for one third of such pairs the longer isoform is less than 38% longer than the shorter one (see Data set 1).

Figure 3 (sim2).

A: DGE detection performance of edgeR applied to three different count matrices (simplesum, scaledTPM, featureCounts), with or without including an offset representing the average transcript length (for simplesum and featureCounts, avetxl indicates that such offsets were used). Including the offset or using the scaledTPM count matrix leads to improved FDR control compared to using simplesum or featureCounts matrices without offset. The curves trace out the observed FDR and TPR for each significance cutoff value. The three circles mark the performance at adjusted p-value cutoffs of 0.01, 0.05 and 0.1. B: stratification of the results in A by the presence of differential isoform usage. The improvement in FDR control seen in A results from an improved treatment of genes with differential isoform usage, while all methods perform similarly for genes without differential isoform usage.

Figure 4. Comparison of log-fold change estimates from edgeR, based on simplesum and scaledTPM count matrices, in four different data sets.

For the simulated data set (sim2), where signals have been exaggerated to pinpoint underlying causes of various observations, genes with induced DTU (whose true overall log-fold change is 0) show a clear overestimation of log-fold changes when using simplesum counts. However, none of the real data sets contain a similar population of genes, suggesting that for many real data sets, simple gene counting leads to overall similar conclusions as accounting for underlying changes in transcript usage.

Discussion

In this article, we have contrasted transcript- and gene-resolution analyses in terms of both abundance estimation and statistical inference, and illustrated that gene-level results are often more accurate, powerful and interpretable than transcript-level results. Not surprisingly, however, accurate transcript-level estimation and inference play an important role in deriving appropriate gene-level results, and it is therefore imperative to continue improving abundance estimation and inference methods applicable to individual transcripts, since misestimation can propagate to the gene level. We have shown that when testing for changes in overall gene expression (DGE), traditional gene counting approaches may lead to an inflated false discovery rate compared to methods aggregating transcript-level TPM values or incorporating correction factors derived from these, for genes where the relative isoform usage differs between the compared conditions. These correction factors can be calculated from the output of transcript abundance programs, using e.g., the provided R package (tximport). It is important to note that the average transcript length offsets must account for the differences in transcript usage between the samples and thus using (sample-independent) exon-union gene lengths will not improve performance.

On the six data sets studied here, simple counting with featureCounts led to very similar conclusions as estimated gene counts from Salmon, when combined with count-based statistical inference tools such as edgeR and DESeq2. Moreover, p-value distributions and mean-variance relationships were similar for actual and estimated counts. Taken together, this suggests that the negative binomial assumption made by the count-based tools is flexible enough to accommodate also estimated counts. All evaluated counting approaches, with and without the inclusion of average transcript length offsets, gave comparable DGE results for genes where DTU was not present. Thus, the extent of the FDR inflation in experimental data depends on the extent of DTU between the compared conditions; notably, our simulation introduced rather extreme levels of DTU, hence the inflated FDR, and the difference between the approaches was considerably smaller in real data sets. Recent studies have also shown that many genes express mainly one, dominant isoform²⁶ and for such genes, we expect that simple gene counting will work well.

All evaluations in this study were performed using well-established count-based differential analysis tools. These methods take as input a matrix of counts, which is assumed to correctly represent the origin of each read in a particular set of libraries. However, due to sequence similarities among transcripts or genes, there is often a hidden uncertainty in the feature abundance estimates, even when the set of input reads is fixed. With the development of fast, alignment-free abundance estimation methods, this uncertainty can now be estimated rapidly using bootstrap approaches (see e.g. Figure 1b). Method development is currently underway in the field to account for this uncertainty in the differential expression analysis (e.g., MetaDiff²⁷, sleuth⁶), which has the potential to improve performance of both DTE and DGE analyses. If such methods are based on (potentially transformed) aggregated transcript counts as gene-level abundance measures, DGE analysis will still be affected by the presence of DTU, and thus could benefit from the inclusion of average transcript length offsets, or by instead using the sum of transcript TPMs as gene abundance measures.

Our results highlight the importance of carefully specifying the question of interest before selecting a statistical approach. Summarization of abundance estimates at the gene level before performing the statistical testing should be the method of choice if the interest is in finding changes in the overall transcriptional output of a gene. However, it is suboptimal if the goal is to identify genes for which at least one of the transcripts show differences in transcriptional output, since it may miss genes where two transcripts change in opposite directions, or where a lowly expressed transcript changes. For gene-level detection of DTE (that is, whether any transcript showed a change in expression between the conditions), statistical testing applied to aggregated gene counts led to reduced power and slightly inflated FDR compared to performing the statistical test on the transcript level and aggregating results within genes (Supplementary Figure 9). Statistical inference on aggregated transcript TPMs (scaledTPM) showed low power for detecting changes that did not affect the overall transcriptional output of the gene, as expected. An alternative to DTE analysis, for potential improved interpretability, is to perform a combination of DGE and DTU analyses, both resulting in gene-level inferences. Table 2 summarizes our results and give suggested workflows for the different types of analyses we have considered.

Finally, we note that abundance estimation at the gene level can reduce the impact of technical biases on expression levels, which have been shown to lead to estimation errors, such as expression being attributed to the wrong isoform²⁸. Non-uniform coverage from amplification bias or from position bias (3’ coverage bias from poly-(A) selection) can result in unidentifiable transcript-level estimation. While correction of technical artifacts in coverage can be attempted computationally, through estimation of sequence- and position-specific biases²⁰, we note that such errors and estimation problems are also minimized when summarizing expression to the gene level. This being said, there may of course be situations where a direct transcript-level analysis is appropriate. For example, in a cancer setting where a specific deleterious splice variant is of interest (e.g., AR-V7 in prostate cancer²⁹), inferences directly at the transcript level may be preferred. However, while this may be preferred for individual known transcripts, transcriptome-wide differential expression analyses may not be warranted, given the associated multiple testing cost.

Data availability

F1000Research: Dataset 1. Data set 1, 10.5256/f1000research.7563.d114722

F1000Research: Dataset 2. Data set 2, 10.5256/f1000research.7563.d114723

F1000Research: Dataset 3. Data set 3, 10.5256/f1000research.7563.d114724

F1000Research: Dataset 4. Data set 4, 10.5256/f1000research.7563.d114725

F1000Research: Dataset 5. Data set 5, 10.5256/f1000research.7563.d114726

F1000Research: Dataset 6. Data set 6, 10.5256/f1000research.7563.d114730

Software availability