Relative Abundance of Transcripts (RATs): Identifying differential isoform abundance from RNA-seq

DTU calling

RATs identifies DTU independently at both the gene and transcript levels using an efficient implementation of the G-test of independence²⁴, without continuity corrections. The criteria RATs uses to identify DTU are described in detail below.

Pre-filtering

Prior to statistical testing by either method, RATs first filters the input isoform abundance data to reduce both the number of low quality calls and the number of tests carried out. Specifically: (i) isoform ratio changes can only be defined for genes that are expressed in both conditions, with at least two isoforms detected, and (ii) transcript abundances must exceed an optional minimum abundance threshold. Transcripts with abundances below the threshold are considered as not detected.

Statistical significance

Significant changes in relative transcript abundance are detected using two separate approaches: one at the gene level and the other at the transcript level. At the gene level, RATs compares the set of each gene’s isoform abundances between the two conditions to identify if the abundance ratios have changed. At the transcript level, RATs compares the abundance of each individual transcript against the pooled abundance of its sibling isoforms to identify changes in the proportion of the gene’s expression attributable to that specific transcript. Both methods include the Benjamini-Hochberg false discovery rate correction for multiple testing²⁵. These tests are performed on the summed abundance of each isoform across the replicates.

Effect size

Transcripts whose absolute difference in isoform proportion is below a set threshold are rejected, even if the difference is statistically significant.

Reproducibility

RATs provides the option to use the bootstrapped abundance estimates obtainable from alignment-free quantification tools to apply a reproducibility constraint on the DTU calls, by randomly selecting individual quantification iterations from each replicate and measuring the fraction of these iterations that result in a positive DTU classification. Typically, each sample is represented by the mean abundance of each transcript, calculated across the quantification iterations. However, this loses the variance information of the quantification. By referring back to the quantification iterations, RATs highlights cases where the quantification was unreliable due to high variability and therefore the DTU result should also be considered unreliable. Similarly, RATs optionally also measures the reproducibility of the DTU results relative to the inter-replicate variation by iteratively sub-setting the samples pool.

Implementation

RATs is implemented in R²⁶ and has been freely distributed through Github as an R source package since August 2016. RATs accepts as input either a set of R tables with abundances (with or without bootstrap information), or a set of Salmon¹⁴ or Kallisto¹¹ output files. An annotation table mapping the correspondence between transcript and gene identifiers is also required, either provided directly or inferred from a GTF file. Results are returned in the form of R data.table objects²⁷. Along with the DTU calls per transcript and gene, the tables record the full provenance of the results. Convenience functions are provided for summary tallies of DTU and isoform-switching results, for ID retrieval, and for visualization of the results via ggplot2 (v2.2.1)²⁸. Details on these are available through the user manual of the package. Once created, all plots produced by RATs remain customisable via standard ggplot2 operations.

Performance

The performance was assessed in two ways. Firstly, the false positives (FP) performance of RATs (v0.6.2) for detection of DTU between two groups relative to the level of experimental replication was measured on groups generated by random selection without replacement from a pool of 16 high-quality wild-type Colombia-0 Arabidopsis thaliana replicates^{29 1}. This was iterated 100 times for each replication level in the range 3 ≤ n ≤ 8. As the two groups are drawn from the same condition, any positive DTU calls must be considered to be false positives. For each iteration, we recorded the fraction of genes and transcripts that were reported as DTU, relative to the total number of genes or transcripts tested in that iteration. The commands and scripts used are from the RATs Github repository.

Secondly, two simulated datasets³⁰ were used to benchmark the sensitivity (s, the fraction of the 1000 DTU events actually detected), false discovery rate (FDR, the fraction of reported DTU events that is not part of the 1000 “real” events) and Matthews correlation coefficient (MCC) of RATs, SUPPA2 and DRIM-Seq. The datasets were made of simulated RNA-seq reads based on the transcriptome annotation and to match realistic RNA-seq transcript expression values. To create the second condition, the abundance values of the two most abundant transcript isoforms originating from a gene locus were swapped for 1000 well-expressed coding gene loci. The transcriptome annotation used for both Human and fly comprised only annotated protein coding genes (13937 in the Drosophila, 20410 in the human) leaving a number of other classifications of gene unaccounted for (1745 in the Drosophila, 41483 in the human). These genes constitute a convenient negative set for simulation and should have no expression, save for any reads misallocated to them by the quantification tools. The simulated datasets were obtained from ArrayExpress² and quantified with both Kallisto (v0.44;¹¹ and Salmon (v0.9.1;¹⁴ using the respective complete annotations that match the simulation of the datasets (Ensembl v70 for the Drosophila and Ensembl v71 for the human;³⁰). The sensitivity, FDR and MCC were measured for a range of comparable parameters between RATs (v0.6.4), SUPPA2 (v2.3) and DRIM-Seq (v1.6, Bioconductor v3.6, R v3.4). No transcript abundance pre-filter was imposed on any of the three DTU tools, and the significance level was set to 0.05 for all runs. For RATs and SUPPA2, three thresholds for the effect size (difference in proportion) were tested; the RATs’ current default of 0.2, and more permissive values 0.1 & 0.05. For DRIM-Seq, threshold values of the likelihood ratio were explored from 0-30. Finally, RATs reproducibility thresholds were explored in the range of 0.8-0.95 for the quantification reproducibility and 0.55-0.85 for the inter-replicate reproducibility. The tool performance was measured using annotations comprised of all annotated genes and only protein coding genes.

Comparison on a real 2-condition dataset

To test the ability of RATs to identify known instances of DTU, we compared it against validated instances of DTU from publicly available RNA-seq data. We took read data from Deng et al. (2013, 31), who identified non-DGE changes in the isoform levels of genes between three human patients with idiopathic pulmonary fibrosis (IPF) and three lung cancer patients used as controls. The dataset contains 25 million 54-base long single-end Illumina reads per lung tissue sample. As in the original at study, we used Ensembl v60³² as the source of the reference human genome and its annotation, in which each of the three discussed genes features two isoforms. Unlike the original study, we used Salmon (v0.7.1, with sequence bias correction enabled, 100 bootstrap iterations and default values for the remaining parameters, using k=21 for the index) to quantify the isoform abundances. DTU was identified by RATs v0.6.2. For comparison, we repeated the quantification and DTU analysis of the data with the same tool versions and parameters, but using the annotation and assembly from Ensembl v87, the current version at the time of this study.

We also submitted the quantification data to SUPPA2, in its psiPerIsoform mode, and to DRIM-Seq. For a fair comparison, we tried to minimize variability in the parameters and data type used by the three tools. As SUPPA2 offered no abundance pre-filtering, RATs and DRIM-Seq were run with abundance threshold values of 0. The p-value cut-off was set at 0.05 for all three tools, using the corrected p-values where available. For the difference in isoform proportion (SUPPA2 and RATs) the threshold was set at 0.20. No threshold was set for the fold-changes in DRIM-Seq. SUPPA2 required and was provided with TPM abundances. For consistency in the use of abundances normalised for transcript length, RATs and DRIM-Seq were also provided with TPM, but the values were scaled up to the average library size of 25M reads, as their testing methods expect counts and would be under-powered if used directly with TPMs. Again, the commands and scripts used are available from the RATs Github repository.

False positives performance

Both the gene-level and transcript-level approaches to identifying DTU implemented in RATs achieved a median FP fraction <0.05 on our A. thaliana dataset, even with only three replicates per condition (Figure 2A). Higher replication results in both a reduction in the number of false positives and restricts the false positives to smaller effect sizes (Figure 2B). The gene-level and transcript-level approaches, however, have different strengths and weaknesses. Simultaneously utilizing the expression information across all the isoforms in a gene makes the gene-level test sensitive to smaller changes in relative expression, compared to testing transcripts individually, but it also makes the gene-level test more prone to false positives. Figure 2 shows that the gene-level test has a higher FP fraction than the transcript-level test, irrespective of replication level or effect size, although the two methods converge for highly replicated experiments or large effect sizes. Furthermore, the gene-level test only identifies the presence of a shift in the ratios of the isoforms belonging to the gene, without identifying which specific isoforms are affected. The transcript-level test, in contrast, directly identifies the specific isoforms whose proportions are changing and has fewer false positives than the gene-level test. However, considering each isoform independently requires a larger number of tests to be performed, thus resulting in a greater multiple testing penalty.

Figure 2. False positives (FP) performance of RATs as a function of replication level.

False positive fraction measured over 100 permutation iterations of randomly selected (without replacement) replicates from a pool of 16 high-quality wild-type Colombia-0 Arabidopsis thaliana replicates from Froussios et al. (2017,²⁹). [A] FP fraction of each bootstrap iteration, for default values of all RATs parameters (v0.6.2), across a range of replication levels, separately for the gene-level test (red) and transcript level test (blue). [B] Mean FP fraction by replication level, as a function of the effect size threshold (effect size = difference between conditions of an isoform’s proportion). For a gene, the effect size is defined as the largest proportion difference observed among that gene’s isoforms. In every iteration, the FP fraction was calculated against the number of genes or transcripts that were eligible for testing each time (a number which remains very stable across iterations and replication levels – see Extended data 1³³).

Comparative performance on simulated DTU

The sensitivity, FDR and MCC performance of RATs, SUPPA2 and DRIM-Seq using Salmon transcript quantifications of annotated protein coding gene isoforms are summarised in Figure 3. Tested with the simulated Human dataset, the parameter defaults for RATs (quantification reproducibility >95%, inter-replicate reproducibility >85% & effect-size >0.2) result in a sensitivity of s = 0.55, MCC = 0.71 and FDR = 0.04, outperforming both other tools. With the same thresholds, SUPPA2 has a higher sensitivity (s = 0.61) but poorer FDR performance (FDR = 0.33). Direct comparison with DRIM-Seq is complicated by different methods for measuring DTU effect-sizes between the tools, however for a likelihood-ratio threshold of 30, DRIM-Seq has similar FDR performance to RATs (FDR = 0.06), but worse sensitivity (s = 0.47). These differences persist for the simulated drosophila dataset. DRIM-Seq consistently shows the lowest sensitivity (≤0.65), while maintaining a FDR ≤0.2 in any of the tried parameter sets. SUPPA2 is the most sensitive of the three tools (0.6 ≤ s ≤ 0.9), but also has the highest FDR (0.35 ≤ FDR ≤ 0.65 in human, 0.10 ≤ FDR ≤ 0.25 in Drosophila). RATs can match the sensitivity of SUPPA2 while maintaining a lower FDR than SUPPA2 by relaxing its quantification reproducibility (Qrep) and inter-replicate reproducibility (Rrep) thresholds. At the highest effect-size thresholds (Dprop_RATs = 0.2 and lr_DRIM-Seq = 0.3) DRIM-Seq has a comparable FDR to that of RATs. Surprisingly, the sensitivity, MCC and FDR of DRIM-Seq is not strongly sensitive to variations in the likelihood ratio effect-size threshold. Consequentially, RATs has worse FDR performance, but better sensitivity than DRIM-Seq at lower effect-size thresholds. Across all the simulated dataset and parameter combinations the gene-level test implemented in RATs shows higher sensitivity and higher FDR compared with the results from the transcript-level test. Extending the test to isoforms from the full set of annotated genes, rather than only those from protein coding genes, adds a considerable number of additional true negatives (Drosophila: 1745, human: 4148, see Section: Performance) resulting in a small increase of FDR and slight reduction of MCC for all tools in both datasets (Extended data 2³³). Similarly, using Kallisto isoform expression quantifications in place of the quantifications from Salmon does not strongly affect the results (Extended data 2³³). The performance results of RATs on these simulated datasets are in good agreement with those presented in Love et al. (2018,³⁴), which also demonstrates that the performance of RATs is similar to, or exceeds, the performace of other DTU tools , including DRIM-seq, SUPPA2 or DEX-Seq.

Figure 3. Performance comparison between RATs, SUPPA2 and DRIM-Seq.

The performance was assessed on the human [A] and Drosophila [B] simulated datasets from ArrayExpress E-MTAB-3766³⁰, over a range of threshold values for the effect size (RATs - Dprop, SUPPA2 - dPSI, DRIM-Seq likelihood ratio - lr) and confidence in the result (RATs quantification reproducibility – Qrep, RATs inter-replicate reproducibility - Rrep). The statistical significance cut-off was at 0.05 for all cases. The measures of performance are the sensitivity, false discovery rate (FDR) and Matthews correlation coefficient (MCC). The datasets were quantified using Salmon 0.9.2 and the metrics were calculated accounting only for the genes strictly listed in the “truth” sets. The results using Kallisto for the quantification are practically identical (see Extended data 2³³).

Recapitulating published validated examples of DTU

After pre-filtering, Deng et al. (2013, 31) tested 3098 Ensembl v60 genes for DTU by quantifying their isoform proportions with RAEM³⁵ and using Pearsons Chi-squared test of independence with a FDR threshold of 5%. They identified 248 genes that were not differentially expressed but displayed significant DTU. Subsequently, they confirmed three of them with qRT-PCR: TOM1L1 (ENSG00000141198), CMTM4 (ENSG00000183723), and PEX11B (ENSG00000131779). Table 1 shows the fraction of the 248 DTU genes identified in this study that were also called by RATs, SUPPA2 and DRIM-Seq, as well as each tool’s verdict on each of the three validated genes. The genes reported as DTU by RATs are listed in Extended data 3 & 4³³ respectively, based on the Ensembl v60 and v87 human annotations.

None of the three tools recapitulated the reported 248 genes well, with the highest fraction of 26% achieved by DRIM-Seq possibly due to a tendency to over-predict (see next section). Of the three validated genes, only CMTM4 is reported by all methods, and only SUPPA2 reports all three genes. Although the rejection of TOM1L1 and PEX11B by DRIM-Seq was due to poor statistical significance, RATs reported that the changes found were both statistically significant and of sufficient effect size. Instead, RATs rejected the genes on the grounds of poor reproducibility (see Section: DTU Calling).

There have been extensive changes in the human transcriptome annotation since Ensembl v60. We hypothesized that these changes could have a significant impact on the set of genes identified in Deng et al. (2013, 31). Table 2 shows that in addition to the new genome assembly, the human transcriptome complexity has increased significantly from Ensembl v60 to the more recent v87. Changing the version of the human annotation from Ensembl v60 to v87 removes 10,253 gene IDs and adds 15,839 new ones. Re-quantifying the RNA-seq data with the updated annotation and re-calling DTU resulted in similarly poor overlap between the tools’ results and the original report (see Extended data 5³³). Of the three validated genes, TOM1L1 was unanimously rejected by all methods, CMTM4 remained unanimously reported as DTU, and PEX11B was reported as DTU by RATs and SUPPA2, but not by DRIM-Seq.

The isoform abundances in Figure 4 reveal that all three genes showed plausible shifts in relative isoform abundance with the Ensembl v60 quantifications, but only PEX11B showed the same shift with Ensembl v87. Instead, TOM1L1 showed no significant changes in any of its 23 isoforms and the primary isoform in the Control samples changed from isoform 2 (ENST00000445275) to isoform 1 (ENST00000348161). CMTM4 shows a similar abundance shift with v87 as it did with v60, but the isoforms implicated changed from isoforms 1 (ENST00000330687) and 2 (ENST00000394106) to isoforms 1 and 5 (ENST00000581487). These changes of context raised questions about the qRT-PCR validation performed in the original analysis of the data³¹. Indeed, when the reported qRT-PCR primers were aligned to the Ensembl v87 sequence and annotation (see Extended data 6³³), only the primers for PEX11B yielded the same conclusion as with Ensembl v60. For TOM1L1, the primers intended for ENST00000445275 no longer matched that isoform, but matched two other isoforms instead (ENST00000570371 and ENST00000575882). Additionally, the primers intended to quantify the gene as a whole failed to match half of the gene’s new isoforms, and the two sets of captured isoforms did not overlap completely and were thus incomparable in any meaningful way. As a consequence, the qRT-PCR intensities measured in the original study are actually impossible to interpret in the context of the updated annotation and the originally reported conclusion is likely wrong. For CMTM4 the primers reported matched multiple but not all isoforms, casting doubt on the interpretation of the qRT-PCR measurements for this gene as well. Only for PEX11B did the primers target the isoforms in a way that would give interpretable results and indeed lead to the same conclusion as originally reported³¹.

Figure 4. Relative abundance of isoforms for the three validated genes from 31, as re-quantified with Salmon 0.7.1¹⁴ using two versions of the Ensembl annotation.

Isoform IDs on the x axis were replaced with simple numbers to minimize clutter, but the mapping of number to ID is maintained between the two annotations. The y axis represents the relative abundance of each isoform. In red are the quantifications from the three replicates of the Control condition, and in blue are those from the IPF condition. The full version of the plots by RATs, including the full isoform IDs, is available in Extended data 5³³.

Comparison of DTU methods against Deng et al. (2013, 31)

Table 3 summarises the results obtained by RATs, SUPPA2 and DRIM-Seq for the Deng et al. (2013, 31) dataset using Ensembl v60 (same as the original study) and Ensembl v87 (current version at time of the present work). With either annotation, DRIM-Seq reported the most DTU genes – almost 1000 with v60 and almost 1700 with v87. The RATs gene-level method reported fewer genes by a factor of 1.5 and 2 respectively compared to DRIM-Seq with each annotation. SUPPA2 reported several hundred transcripts more than RATs, but at the gene level the numbers were comparable. RATs and DRIM-Seq reported more genes and transcripts with v87 of the annotation than with v60, whereas SUPPA2 reported slightly fewer with v87. Despite overall similar volume of results between the two versions of the annotation, it is evident from Table 3 that the overlap of the results between annotations is poor for all methods. For RATs and SUPPA2, only 30–40% of the genes reported with Ensembl v60 were also reported with v87. For DRIM-Seq this overlap was 55% of its Ensembl v60 results.

The overlap of results between different methods is similar to the overlap of results between annotations, as shown in Table 4. 97% of the genes reported by gene-level method in RATs are also identified as DTU by the transcript-level method. Among all the pairwise comparisons of RATs, SUPPA2 and DRIM-Seq, however, the highest level of agreement at both transcript and gene level is between SUPPA2 and RATs. SUPPA2 identifies DTU in 53% of the transcripts that are called as DTU by the transcript-level method in RATs, however RATs calls DTU for only 35% of the transcripts identified as DTU by SUPPA2. DRIM-Seq consistently reports a higher number of DTU identifications than either RATs or SUPPA2, but still only manages at most 43% agreement with the other two tools.

RATs and SUPPA2 are more similar than implied by the level of agreement presented in Table 4. Figure 5 shows that the novel reproducibility testing feature in RATs, which discounts DTU identification from highly variable quantifications (see Section: DTU Calling), is responsible for rejecting 43% of the SUPPA2 DTU transcripts and 28% of the DRIM-Seq DTU genes that pass the significance and effect size filtering criteria. 53% of the DRIM-Seq results and, perplexingly, 18% of the SUPPA2 results are rejected due to the effect size filter (after passing the significance testing, but prior to the reproducibility filter), despite all the tools operating on the same input isoform quantifications.

Figure 5. Causes of rejection by RATs of results reported by SUPPA2 or DRIM-Seq, expressed as proportion of the total DTU identifications reported by SUPPA2 (1252 transcripts) or DRIM-Seq (1680 genes).

The colours represent the different criteria imposed by RATs. Since no abundance pre-filtering was enabled for any of the tools, there are no rejections caused by the transcript abundance and the effective number of expressed isoforms. DTU: Differential transcript usage.

Hardware requirements and run times

RATs’ runtime and memory consumption depend on the size of the annotation and the number of bootstraps iterations. Where multiple processing cores are available, RATs can be instructed to take advantage of them. The runtime and maximum memory usage for the two simulated datasets from our benchmarks, running on a high-specification laptop, are shown in Table 5.

Underlying data

The Arabidopsis thaliana RNA-sequencing data used in this study is available from ArrayExpress under the study E-MTAB-5446. The simulated Homo sapiens and Drosophila Melanogaster datasets are available from ArrayExpress under the study E-MTAB-3766. The Deng et al. (2013, 31) data are available from the European Nucleotide Archive, or the Short Read Archive, under the study SRA048904.

Extended data

Extended data are available along with the source code from GitHub and archived with Zenodo

Zenodo: Extended data. bartongroup/RATS: RATs 0.6.5 - R source package, http://doi.org/10.5281/zenodo.2556564³³

Licence: MIT

Software availability

The RATs R package is open source and available through Github

Source code: https://github.com/bartongroup/RATs.

Archived source code: http://doi.org/10.5281/zenodo.2556564³³

Licence: MIT