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

The biological importance of changes in RNA expression is reflected by the wide variety of tools available to characterise these changes from RNA-seq data. Several tools exist for detecting differential transcript isoform usage (DTU) from aligned or assembled RNA-seq data, but few exist for DTU detection from alignment-free RNA-seq quantifications. We present the RATs, an R package that identifies DTU transcriptome-wide directly from transcript abundance estimates. RATs is unique in applying bootstrapping to estimate the reliability of detected DTU events and shows good performance at all replication levels (median false positive fraction < 0.05). We compare RATs to two existing DTU tools, DRIM-Seq & SUPPA2, using two publicly available simulated RNA-seq datasets and a published human RNA-seq dataset, in which 248 genes have been previously identified as displaying significant DTU. RATs with default threshold values on the simulated Human data has a sensitivity of 0.55, a Matthews correlation coefficient of 0.71 and a false discovery rate (FDR) of 0.04, outperforming both other tools. Applying the same thresholds for SUPPA2 results in a higher sensitivity (0.61) but poorer FDR performance (0.33). RATs and DRIM-seq use different methods for measuring DTU effect-sizes complicating the comparison of results between these tools, however, for a likelihood-ratio threshold of 30, DRIM-Seq has similar FDR performance to RATs (0.06), but worse sensitivity (0.47). These differences persist for the simulated drosophila dataset. On the published human RNA-seq dataset the greatest agreement between the tools tested is 53%, observed between RATs and SUPPA2. The bootstrapping quality filter in RATs is responsible for removing the majority of DTU events called by SUPPA2 that are not reported by RATs. All methods, including the previously published qRT-PCR of three of the 248 detected DTU events, were found to be sensitive to annotation differences between Ensembl v60 and v87.

High-throughput gene regulation studies have focused primarily on quantifying gene expression and calculating differential gene expression (DGE) between samples in different groups, conditions, treatments, or time-points. However, in higher eukaryotes, alternative splicing of multi-exon genes and/or alternative transcript start and end sites leads to multiple transcript isoforms originating from each gene. Since transcripts represent the executive form of genetic information, analysis of differential transcript expression (DTE) is preferable to DGE. Unfortunately, isoform-level transcriptome analysis is more complex and expensive since, in order to achieve similar statistical power in a DTE study, higher sequencing depth is required to compensate for the expression of each gene being split among its component isoforms. In addition, isoforms of a gene share high sequence similarity and this complicates the attribution of reads among them. Despite these challenges, several studies have shown that isoforms have distinct functions [1][2][3] and that shifts in individual isoform expression represent a real level of gene regulation 4-7 , suggesting there is little justification for choosing DGE over DTE in the study of complex transcriptomes.
It is possible to find significant DTE among the isoforms of a gene, even when the gene shows no significant DGE. This introduces the concept of differential transcript usage (DTU), where the abundances of individual isoforms of a gene can change relative to one another, with the most pronounced examples resulting in a change of the dominant isoform (isoform switching). The definitions of DGE, DTE and DTU are illustrated in Figure 1.
To quantify the isoforms and assess changes in their abundance, most existing tools for DTE and DTU analysis (e.g. Cufflinks 8 , DEXSeq 9 , LeafCutter 10 ) rely on reads that either span splicejunctions or align to unique exons. However, with the newest generation of transcript quantification tools (Kallisto 11,12 , Sailfish 13 , Salmon 14 ), reads are aligned to neither the transcriptome nor the genome. Instead, these tools combine a pseudo-mapping of the k-mers present within each read to the k-mer distributions from the transcriptome annotation with an expectation maximization algorithm, to infer the expression of each transcript model directly. Such alignment-free methods are much faster than the traditional alignment-based methods (RSEM 15 , TopHat2 16 , STAR 17 ) or assembly-based methods (Cufflinks 8 , Trinity 18 ), making it feasible to repeat the process many times on iterative subsets of the read data and, thus, quantify the technical variance in the transcript abundance estimates. However, the lack of alignments prevents these new methods from being compatible with differential expression methods such as Cufflinks, DEXSeq and Leafcutter. Instead, Sleuth 19 is a tool that handles DTE analysis from alignment-free transcript quantifications. DTU analysis is currently less straight-forward. SwitchSeq 20 focuses on a particular subset of DTU analysis from alignment-free data, namely isoform switching, whereas iso-kTSP 6 identifies both DTU and isoform switching, but focuses on the highest-ranking pair of change-exhibiting isoforms per gene. SUPPA 21,22 , on the other hand, primarily The expression of two genes (Gene A and Gene B), with 3 and 2 isoforms respectively, is compared across two conditions (Condition 1 and Condition 2). The horizontal width of each coloured box represents the abundance of the relevant gene or transcript. A negative differential expression result (red cross-mark) for a given entity in any one of the three analysis types does not exclude that same entity from having a positive result (green tick-mark) in one of the other two analysis types. The relative isoform abundances in [iii] are scaled to the absolute isoform abundances in [ii], which in turn are scaled to the gene expressions in [i]. Gene A is differentially expressed, but only two of its three isoforms are differentially expressed (A.2 and A.3). Proportionally, Gene A's primary isoform (A.3) remains the same, but the ratios of the two less abundant isoforms change. Gene B is not differentially expressed, but both its isoforms are differentially expressed, and demonstrate an example of isoform switching. DGE: Differential gene expression, DTE: Differential transcript expression, DTU: Differential transcript usage.
identifies differential splicing events at the junction level, with recent developmental versions having added isoform-level capability. Finally, DRIM-Seq 23 identifies DTU directly from quantification data, but defines the effect size as a fold change which may not be the most appropriate way to compare proportions.
In this paper, we present RATs (Relative Abundance of Transcripts), an R package for identifying DTU directly from isoform quantifications. It is designed to use alignment-free abundance data and is the only tool that exploits bootstrapping to assess the robustness of the DTU calls. RATs provides raw, summary and graphical results, allowing for ease of use as well as for advanced custom queries, and the R language is the environment of choice for many widely-used DGE and DTE tools, allowing for easy integration of RATs in existing workflows. We assess the accuracy of RATs in comparison to SUPPA2 and DRIM-Seq and find RATs to perform at as well as or better than its competitors. Finally, we demonstrate that the results of both RNA-seq based and qRT-PCR based analyses are sensitive to the annotation used for transcript quantification and primer design, respectively.

DTU calling
RATs identifies DTU independently at both the gene and transcript levels using an efficient implementation of the G-test of independence 24 , 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 25 . 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 26 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 14 or Kallisto 11 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.

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 30 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

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 singleend Illumina reads per lung tissue sample. As in the original at study, we used Ensembl v60 32 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 genelevel 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.
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%, interreplicate 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 effectsizes 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 2 https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3766/ 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 33 ). 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 33 ).
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 33 ). Similarly, using Kallisto isoform expression quantifications in place of the quantifications from Salmon does not strongly affect the results (Extended data 2 33 ). The performance results of RATs on these simulated datasets are in good agreement with those presented in Love et al. (2018,34 ), 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.

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 35 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 33 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 33 ). 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 31 . Indeed, when the reported qRT-PCR primers were aligned to the Ensembl v87 sequence and annotation (see Extended data 6 33 ), 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 31 . 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.

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

SUPPA2 SUPPA2
(tr. aggr.)   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.

Discussion
Reliable identification of differential isoform usage depends critically on i) the accuracy of the upstream isoform expression quantifications, and ii) on the accuracy of the annotation they use. RATs is the first differential isoform usage tool to include the reproducibility of the upstream isoform expression quantifications to refine its DTU identifications, directly addressing the accuracy of the upstream isoform expression quantifications. Leveraging the bootstrapped isoform expression quantifications from fast modern alignment-free isoform expression quantification tools (such as Kallisto and Salmon) allows RATs to reject those cases of DTU that are based on highly uncertain isoform quantifications. Existing tools rely on the mean isoform abundances, which can hide a large degree of variability, and are thus insensitive to this reproducibility criterion. We recommend running RATs, and the underlying alignment-free isoform expression quantification tools that generate the data it operates on, with at least 100 bootstrap iterations. In the future, experiment-specific transcriptome annotations could be obtained by including a parallel set of full-length isoform RNA-seq data in the experimental design, such as via PacBio sequencing or Oxford Nanopore Direct RNA-seq. An advantage of this approach is that it would better define the transcriptome for the specific experiment 38-41 . This may be of importance for experiments focusing on specific tissues or developmental stages of an organism, where the active transcriptome is likely to be only a subset of the global reference transcriptome of the organism. The authors present a new method for detection of DTU from RNA-seq data, which uniquely leverages quantification uncertainty in the form of inferential replicates. I am not aware of other methods specifically designed to detect DTU as opposed to change in total expression level of the gene, which take into account quantification uncertainty. It is therefore a useful contribution to the methods literature. The authors have taken some length to assess their method against other popular methods on real and simulated datasets, and investigating individual genes with qRT-PCR validation in detail.

Data availability
I have some concerns about the conclusions from the evidence provided in the article, and additionally have requests for further details about the methods, which should be presented in the article itself.

Major comments:
The methods are not sufficiently described, I have the following questions:

1)
What is the input to RATs? Is it TPM or counts or scaledTPM? Should the library size differences be removed prior to providing to RATs or does RATs take care of library size differences internally?
Can the methods described all analyze the same type of experiment, are they all restricted to two-group analyses? Can any of them control for batch effects?
What are the default pre-filtering and post-filtering settings? What is the default minimum abundance threshold or proportion threshold for an isoform to be considered expressed? What is the default effect size cutoff, and how is it implemented per isoform, per gene? What is the default fraction for determining that evidence of DTU is not substantiated across inferential replicates? Likewise, what default fraction for biological replicate variation? I didn't understand why abundance thresholds were not used, as described here, "

2)
No transcript ," and also " abundance pre-filter was imposed on any of the three DTU tools As SUPPA2 offered no " abundance pre-filtering, RATs and DRIM-Seq were run with abundance threshold values of 0.

As shown in Soneson
(2016 ) and Love (2018 ), performance of a number of DTU methods is et al.
et al. greatly improved by filtering out lowly expressed transcripts. It can be inferred from the title of the former paper: "Isoform prefiltering improves performance of count-based methods for analysis of differential ". transcript usage 1 2 ". transcript usage SUPPA2 does have an abundance pre-filtering option, which was used in Love (2018 ): et al. "We enabled a filter to remove transcripts with less than 1 TPM. TPM filtering is a command-line option available during the diffSplice step of SUPPA2 and this greatly improved the running time without loss of ". sensitivity From the SUPPA2 manual: "-th | --tpm-threshold: Minimum expression (calculated as average TPM value " within-replicates and between-conditions) to be included in the analysis. (Default: 0).
Given that all the methods have abundance and/or proportion filters available, that filters are recommended by at least two of the three methods in their documentation (DRIMSeq and RATs), and at least two independent review papers (not introducing methods) have shown that abundance and/or proportion filtering improves performance of methods, I can't see why the choice was made to not use filters.
It is mentioned in the false positive analysis that the median FP fraction was less than 0.05 and a 3) horizontal line is drawn on Fig 2A and B. This is misleading, as the adjusted p-values are being thresholded at 0.05 (I assume), and in a null comparison the rate of false positives from an adjusted p-values should be 0, not 0.05. Drawing or mentioning a 0.05 cutoff would be relevant for the p-values (uncorrected), but has no bearing on the adjusted p-values. This may confuse readers.
The authors repeatedly refer to the reported effect size in DRIMSeq being an issue for comparison 4) across methods, e.g. "Direct comparison with DRIM-Seq is complicated by different methods for ", but this is only an issue to the extent that the authors wish measuring DTU effect-sizes between the tools to perform post-hoc filtering on effect size. It is not an issue for null hypothesis testing without post-hoc filtering, because all methods are testing against the null that the underlying proportion of expression across isoforms has the same distribution for control samples and treated samples. However, I agree that for post-hoc filtering, one may want to filter the methods in a similar manner. It should be easy to filter the DRIMSeq results directly on absolute difference in isoform proportion, for example in Love (2018 ) et al. we performed post-hoc filtering for DRIMSeq on the SD of proportions across all samples using a 6-line R function.
As the likelihood ratio statistic should be 1-1 and monotonic with the p-value for DRIMSeq (if the degrees of freedom is constant across genes or transcripts), then I would not compare effect size filtering with likelihood ratio filtering, as the latter is simply filtering the p-value at a lower threshold.

Minor comments:
In the Introduction, the authors state "there is little justification for choosing DGE over DTE in the ". The authors imply that gene-level and transcript-level analysis study of complex transcriptomes are mutually exclusive analyses, when they are not, and so I would suggest to reword or reconsider this statement. I and others have encouraged assessment of total changes in gene expression (DGE) as well as changes in isoform proportion (DTU), as both may be present in an experiment and both may be of biological importance to the system being studied. DGE has the property that the majority of inferential uncertainty which exists in an RNA-seq sample is removed (because it occurs across isoforms within genes), leaving inferential uncertainty from reads mapping across gene loci, but this property of reduced uncertainty does not preclude a transcript-level analysis. 2 2 gene loci, but this property of reduced uncertainty does not preclude a transcript-level analysis. While DTE has advantages, the above sentence claiming that DGE has none overstates a more complex situation in my opinion.
Throughout the paper, the authors refer to "DRIMSeq" as "DRIM-Seq" which is minor but different than the software and publication.
For what it's worth, the transcript-level test is similar conceptually to the current implementation of testForDEU() in DEXSeq which compares the expression of each feature to the sum of expression from all other features of the gene (this is also different from the test described in the original DEXSeq publication). Running DEXSeq on transcript estimated counts with testForDEU() was tested on simulated data in Soneson It is stated that, "The performance results of RATs on these simulated datasets are in good ". However, this seems to be not clearly the agreement with those presented in Love et al (2018) case, which may be due to differences in the simulated data in the two articles, or some other reason. In the present article, DRIMSeq is reported as having lower sensitivity with lower achieved FDR than other methods, SUPPA2 has higher sensitivity and higher FDR, and RATs with various filter thresholds falls in between. In Love (2018 ), DRIMSeq had the opposite performance: et al. higher sensitivity but higher FDR relative to SUPPA2 and RATs run with default filters. However interpretation is made difficult by all the filtering options in Figure 3. It would be easier to compare perhaps if an additional supplementary plot to Figure 3 was made with only the default filter thresholds instead of the filter threshold ranges for all methods. The main commonality across the two benchmarks seems to be that RATs can achieve higher sensitivity than SUPPA2 while maintaining the same precision, for the 5% nominal FDR threshold.
This sentence needs to be made more specific, or else it could be misleading: "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. Specifically which conclusion is likely wrong? From the analysis, it seemed like there is not a problem with the original qRT-PCR intensities and interpretation for at least one of the three genes.
Why is it perplexing that " ". I 18% of the SUPPA2 results are rejected due to the effect size filter didn't follow the authors in that statement.
It is stated, " ...". This implies that the mean of Existing tools rely on the mean isoform abundances inferential replicates is used for statistical testing. It's perhaps subtly different, other methods are typically using the maximum likelihood estimate, which may be different than the mean of the bootstrap distribution, and different than the mean of the Gibbs sampling distribution. I would just say that other tools do not make use of inferential replication. 'RATs' is addressing an important problem: how to quantify changes in transcript isoform usage. Other tools, like 'Sleuth', address a related problem, which is differential transcript expression. Both RATs and Sleuth take advantage of the bootstrapping data that tools like 'kallisto' and 'Salmon' generate when quantifying transcript isoform abundance. By taking advantage of bootstrapping, such tools can estimate the technical variation within the data, to better look for differential changes. Transcript isoform usage is often linked to changes in alternative splicing or isoform specific decay rates (e.g. from NMD). Therefore, having a tool to accurately find changes in isoform usage is vital to our ability to address a range of biological problems.
I have tried a version of RATs. I found that it was easy to install and easy to use. Being able to install bioinformatics software is no guarantee (Mangul , 2018 ). A bonus of RATs is that several figures et al.
can be generated from the data within the tool. This was simple to do, but allowed for you to visualize your 1 can be generated from the data within the tool. This was simple to do, but allowed for you to visualize your data in a straightforward but powerful way. This is one of the rare tools that just works and was relatively intuitive and well-documented.
The paper uses sensible approaches to compare RATs to other tools, including 'SUPPA' and 'DRIMSeq' and the authors found that RATs performed at a similar or better level than the other tools.
One minor point that could be better explained is how RATs uses the bootstrapping data. Does it use it to simply throw out highly variable genes (decreasing FP rate) or does it help get closer to the true rate of biological variation, thus increasing the true positive rate?
In the methods, it would be good to see more explanation on how RATs does its pre-filtering. For example, if a transcript has zero expression in one treatment but a modest to high expression in the other treatment, would RATs keep this transcript or discard it? This would be of interest to people working on RNA decay pathways, such as NMD.
One thing that I would love to see is a comparison of RATs to the DEXSeq/DRIMSeq approach used to address differentiation transcript isoform usage (Love , 2018 ). This tool appeared to also perform et al. well in the publication of this approach, where they used sim data. Therefore I think a comparison to RATs here, using sim and real data (human) would be appropriate.
Reviewer Expertise: I have research experience in wet lab biology and dry lab (computational analysis). I have used many computational tools like RATS for analyses with real data.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
28 February 2019 Referee Report https://doi.org/10.5256/f1000research.19594.r44979

Sophie Shaw
Centre for Genome Enabled Biology and Medicine, University of Aberdeen, Aberdeen, UK

Froussios
have presented here a new tool, RATS, for the identification of differential transcript usage et al. from transcript abundance estimates. RATs was benchmarked and compared to the existing tools DRIM-Seq and SUPPA2 across four different datasets. False positive rate, false negative rate, sensitivity and Matthews correlation coefficient were all measured. When considered as a whole, RATs was found to outperform the other tools. Differing results due to the version of the reference genome used are also discussed. This is a nicely presented manuscript, with well thought out comparisons. The tool will make a good addition to existing RNA sequencing analysis pipelines, especially as the field moves towards alignment free methods.
The rationale for the development of this tool is clearly stated, as there are only a few tools which carry out DTU detection from alignment-free RNA-seq quantifications. The majority of existing tools for DTE and DTU are designed for use with alignment-and assembly-based methods. Of the existing tools described, each has specific uses, and RATs has been presented as a broad "differential transcript usage" identification tool.
The methods of the analysis have been described well, and overall are technically sound. I would like to see an expansion on the description of the statistical method underpinning RATs. Although G-test of independence is cited, a brief description of what this entails and how it differs from existing tools would aid in the understanding of how the tool functions.
However, I have some suggestions concerning the comparison of tools and the datasets selected. With regards to the selection of tools for comparison, SwitchSeq and iso-KTSP are mentioned within the introduction as being able to use transcript abundance estimates, however are not compared to. I assume that this is because they are too specialist in their identification of differential transcript usage and/or isoform switching, but I think the decision to not compare to these tools should be more explicit. The authors have not mentioned the recent pre-print from Cmero et al. (2019) which discusses the development of methods for DTU detection from alignment-free datasets using equivalence classes. The paper uses the same simulated datasets for benchmarking of the tool, and should be considered as another tool to compare to RATs. If this is not deemed as an equivalent method, it should at least be discussed in this manuscript.
With regards to the datasets tested, the published human data set which is shown here is not directly confirming the accuracy of RATs, as the authors show that the qPCR validation within the original study may be inaccurate, and underlying issues are present due to the reference genome version. Although the dataset is being used to compare RATs to SUPPA and DRIM-Seq, it is not validating the tool. I think that this manuscript would benefit from comparison of the three tools using another "real-life" dataset, which has been validated in some way, to support that RATs is detecting known DTU. 1 I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias You can publish traditional articles, null/negative results, case reports, data notes and more The peer review process is transparent and collaborative Your article is indexed in PubMed after passing peer review Dedicated customer support at every stage For pre-submission enquiries, contact research@f1000.com