Evaluation and comparison of computational tools for RNA-seq isoform quantification

Abstract

Background: Alternatively spliced transcript isoforms are commonly observed in higher eukaryotes. The expression levels of these isoforms are key for understanding normal functions in healthy tissues and the progression of disease states. However, accurate quantification of expression at the transcript level is limited with current RNA-seq technologies because of, for example, limited read length and the cost of deep sequencing.

Results: A large number of tools have been developed to tackle this problem, and we performed a comprehensive evaluation of these tools using both experimental and simulated RNA-seq datasets. We found that recently developed alignment-free tools are both fast and accurate. The accuracy of all methods was mainly influenced by the complexity of gene structures and caution must be taken when interpreting quantification results for short transcripts. Using TP53 gene simulation, we discovered that both sequencing depth and the relative abundance of different isoforms affect quantification accuracy CONCLUSIONS: Our comprehensive evaluation helps data analysts to make informed choice when selecting computational tools for isoform quantification.

Keywords: Data analysis; Isoform; Kallisto; Quantification; RNA-seq; RSEM; Salfish; Salmon.

Fig. 1

Workflow for transcript isoform quantification. Sequencing reads were either mapped by STAR aligner or directly fed into alignment-free methods, Salmon, Sailfish or Kallisto. The transcriptome BAM files were quantified by Salmon_aln, eXpress, RSEM or TIGAR2. The genome BAM files were quantified by Cuffquant and then Cuffnorm from the Cufflinks package. The results are summarized into counts and TPM tables for comparison

Fig. 2

Comparisons of the overall performance among different methods and the impact of the number of transcripts on the accuracy of isoform quantification. a Pearson correlation coefficient. b mean absolute relative differences and c-d) The above metrics were broken into separate groups according to the number of annotated transcript isoforms for each gene. The number of transcripts in each group is shown in figure legends. The accuracy metrics were calculated by comparing the estimated counts with the “ground truths” in simulated dataset

Fig. 3

Inconsistency in effective length calculation among methods for short transcripts. a Read counts were estimated correctly for the transcript SNHG25–002. b Methods showed disagreement in estimating TPM values for the same transcript. c The relationship between the effective transcript length estimated by each method and the corresponding transcript length. Only transcripts with length less than 400 nt are shown. Note the transcript length at x-axis is the total number of nucleotides of the transcript in the Gencode Release v25

Fig. 4

Correlation of estimated TPM values for all transcripts between technical replicates of experimental datasets. a UHRR-C1 (x-axis) and UHRR-C2 (y-axis). b HBRR-C4 (x-axis) and HBRR-C6 (y-axis). The R ² value is shown in each figure. Note, x and y-axes represent log₂ transformed estimated transcript TPM values

Fig. 5

Pairwise correlation of estimated TPM values for all transcripts between methods for the HBRR-C4 sample. The distribution of transcripts’ TPMs from each method was plotted on the diagonal panels. Pairwise density plots and R ² values are shown in the lower and upper triangular panels, respectively. R ² values over 0.9 are in bold. Methods are grouped using hierarchical clustering

Fig. 6

Significant difference in estimated read counts for transcript RPS28P7–001 resulting from STAR aligner. A total of 154 reads for RPS28P7–001 were simulated. a The estimated read counts from all eight methods are shown, and they are severely underestimated by the methods using STAR aligner. b The read coverage profiles (coloured in red) in RPS28P7–001 and RPS28–001. The peak paired-end read counts (both ends counted) are shown in brackets. Only a small fraction of reads were mapped back to the RPS28P7 region while the majority of reads were incorrectly mapped to the RPS28 gene

Fig. 7

The impact of sequencing depth and relative abundance on the accuracy of isoform quantification. a Structures of six canonical transcripts of the TP53 gene, and their corresponding identifier in GENCODE v25. b The accuracy of isoform quantification with each of the seven methods under each simulation condition. MARDS was calculated using known and estimated read counts