Source-linked AI summary
Sailfish: Alignment-free Isoform Quantification from RNA-seq Reads using Lightweight Algorithms
Rob Patro, Stephen M. Mount, Carl Kingsford
TL;DR
RNA-seq analysis requires faster ways to quantify expression accurately. Sailfish uses lightweight indexing and counting structures, and remains accurate while running over an order of magnitude faster than existing alternatives.
Problem
Existing tools were difficult to compare because they expose many parameter options, motivating evaluation of Sailfish against alternative expression-quantification tools.
Method
Sailfish uses a lightweight approach based on indexed reference sequences and counting structures for quantification.
Results
Over an order of magnitude faster than existing alternatives, Sailfish retained accurate expression estimates even when competitors sacrificed accuracy for speed.
Takeaways & Limitations
Sailfish’s lightweight indexing and counting structures require substantially less storage than alignment files, totaling 2.4Gb in one quantification setup.
Takeaways & Limitations
Sailfish had not implemented or tested bias correction for some features.
Abstract
from arXiv · showhide
RNA-seq has rapidly become the de facto technique to measure gene expression. However, the time required for analysis has not kept up with the pace of data generation. Here we introduce Sailfish, a novel computational method for quantifying the abundance of previously annotated RNA isoforms from RNA-seq data. Sailfish entirely avoids mapping reads, which is a time-consuming step in all current methods. Sailfish provides quantification estimates much faster than existing approaches (typically 20-times faster) without loss of accuracy.
Methods
Sailfish quantifies transcript abundance without read mapping by indexing reference-transcript k-mers, counting observed read k-mers efficiently, and estimating abundances with an accelerated EM procedure. It also applies optional post-estimation bias correction using transcript-level sequence features.
- Indexing: Sailfish builds a k-mer index containing a BDZ minimum perfect hash function, k-mer counts, and forward and reverse transcript lookup tables.The index depends only on the reference transcripts and k-mer length, so it must be recomputed only when either changes.
- Quantification: Sailfish counts reference k-mers observed in RNA-seq reads using a precomputed perfect hash function and a compact count array.It checks forward and reverse-complement k-mers and uses a majority-rule heuristic to select which orientation contributes counts.
- Quantification: Sailfish groups k-mers into equivalence classes defined by shared transcript membership and frequency, then allocates class counts rather than full k-mer counts.This reduces the representation used during abundance estimation while preserving the transcript relationships encoded by the classes.
- Estimating abundances via an EM algorithm: The abundance estimator alternates E-steps allocating observed equivalence-class counts among transcripts and M-steps updating relative transcript abundances.Transcript length is adjusted as l̂_i = l_i − k + 1 in the M-step.
- Estimating abundances via an EM algorithm: Sailfish replaces standard EM updates with SQUAREM, which uses three successive EM steps to estimate update direction and step size, substantially improving convergence speed.The procedure makes relatively large updates and is repeated for 30 steps in all reported experiments.
- Bias Correction: Sailfish performs bias correction after abundance estimation using random-forest regression on transcript length, GC content, and dinucleotide frequencies.The correction subtracts predicted bias contributions from original RPKMs and can be disabled with --no-bias-correction.
Sailfish: Alignment-free Isoform Quantification from RNA-seq Reads using · Supplementary Figure 1: Effect of k-mer length on retained data and k-mer
Sailfish is introduced as an alignment-free method for quantifying annotated RNA isoforms from RNA-seq data. In the synthetic dataset, increasing k-mer length from 15 to 25 slightly reduced data retention while affecting the proportion of unique k-mers.
- Supplementary Figure 1: Effect of k-mer length on retained data and k-mer: 15–25 was the tested k-mer-length range for processing the synthetic dataset.
- Supplementary Figure 1: Effect of k-mer length on retained data and k-mer: Longer k-mer lengths produced a slight decrease in data retention.Data retention was defined as the ratio of hashable read-set k-mers to all k-mers appearing in the reads.
- Supplementary Figure 1: Effect of k-mer length on retained data and k-mer: Data retention was measured by the ratio of hashable k-mers from the read set to the total number of k-mers in the reads.
- Supplementary Figure 1: Effect of k-mer length on retained data and k-mer: The ratio of unique k-mers in the transcript set was examined simultaneously as k-mer length varied.Unique k-mers were defined as k-mers having a unique locus of origin.
Supplementary Figure 2: Speed of counting indexed k-mers
Supplementary Figure 2 measures k-mer counting and transcript-abundance quantification speed on an 81M-read dataset as concurrent hashing threads increase. Processing time decreases from 34 minutes 26 seconds with one thread to 1 minute 28 seconds with 32 threads.
- Supplementary Figure 2: Speed of counting indexed k-mers: 1 minute and 28 seconds suffices to count all k-mers and quantify transcript abundance in the 81M-read SRX016366 dataset using 32 processing threads.With a single thread, the same processing takes 34 minutes and 26 seconds.
Supplementary Note 1: Additional benefits of the Sailfish approach
Sailfish requires substantially smaller indexing and counting structures than the indexing and alignment files used by most other methods. Its index files depend on unique k-mers and transcriptome k-mer complexity rather than read count.
- Scaling behavior: Sailfish index files grow with the number of unique k-mers and transcriptome k-mer complexity, independent of the number of reads.This contrasts with alignment files, which grow with the number of reads.
Supplementary Figure 3: Correlation plots with qPCR on human brain · tissue and synthetic data
Supplementary Figure 3 compares RSEM, eXpress, and Cufflinks using correlation plots against qPCR-based expression estimates for human brain tissue and true transcript abundance in synthetic data.
- Supplementary Figure 3: Correlation plots with qPCR on human brain: The figure presents correlation plots for RSEM, eXpress, and Cufflinks using the data from Fig. 2.Each plot column is labeled by the method whose output generated that column.
- Supplementary Figure 3: Correlation plots with qPCR on human brain: Each column corresponds to one quantification method’s output.The columns are labeled with the method used to produce their plots.
- Supplementary Figure 3: Correlation plots with qPCR on human brain: The top row evaluates computed RPKM against qPCR-based expression estimates for human brain tissue.These plots assess agreement between the computed RPKM values and qPCR-based estimates.
- tissue and synthetic data: The bottom row evaluates computed RPKM against each transcript’s true abundance in the synthetic dataset.These plots compare computed RPKM with the known transcript abundances in synthetic data.
- Supplementary Figure 3: Correlation plots with qPCR on human brain: Human brain tissue correlations use qPCR-based expression estimates as the reference.The tissue analysis is shown in the figure’s top row.
- tissue and synthetic data: Synthetic-data correlations use true transcript abundance as the reference.The synthetic-data analysis is shown in the figure’s bottom row.
Supplementary Figure 4: Correlation with qPCR on universal human · reference tissue
Supplementary Figure 4 evaluates four methods on a second MACQ dataset from mixed Universal Human Reference tissues, using approximately 93 million reads and the same reference transcripts as the main-text comparison. The methods’ relative accuracy and performance were similar to those observed in the other MACQ dataset.
- reference tissue: Supplementary Figure 4 evaluates four methods on a second dataset from the MACQ study using reads from SRA accession SRX016367.The dataset contains approximately 93 million reads from a mixture of different tissues, termed the Universal Human Reference (UHR).
- reference tissue: The same set of reference transcripts was used as in Figure 2 of the main text.
- reference tissue: The methods’ relative accuracy and performance were similar to those observed in the other MACQ dataset.The passage identifies Sailfish, eXpress, and Cufflinks before truncation.
Supplementary Note 2: Additional details of accuracy analysis
Accuracy was assessed with correlation, RMSE, and median percentage error, which capture complementary aspects of agreement between predicted and true abundances. Low but non-zero outlier estimates from eXpress and Cufflinks degraded Pearson correlation and were filtered using a 0.01 RPKM cutoff.
- Metrics: Accuracy was evaluated using Pearson and Spearman correlations, RMSE, and median percentage error.These metrics gauge accuracy from different perspectives.
- Metrics: Pearson correlation measures captured trends, whereas log-scale correlation discounts zero or very low abundances and RMSE includes zero true or estimated abundances.The metrics therefore differ in which transcripts contribute to the assessment.
- Outlier filtering: eXpress and Cufflinks produced a few very low but non-zero abundance outliers that significantly degraded Pearson correlation.These outliers were discarded by setting estimated RPKM values less than or equal to 0.01 to zero; the cutoff removed outliers without apparently discarding other truly expressed transcripts.
Supplementary Note 3: Parameters for simulated data
Simulated RNA-seq data were generated with FluxSimulator v1.2 using specified expression, fragmentation, reverse-transcription, amplification, and sequencing parameters. The setup used 5,000,000 molecules and 150,000,000 paired-end 76-base reads.
- Simulation framework: FluxSimulator v1.2 generated the simulated RNA-seq data using hg19 annotations and genomic sequences.The reference file was hg19_annotations.gtf, with GEN_DIR set to hg19/chrs.
- Expression: Expression simulation used NB_MOLECULES 5000000 and excluded noncoding transcripts with LOAD_NONCODING NO.TSS_MEAN was 50, while POLYA_SCALE and POLYA_SHAPE were NaN.
- Sequencing and library preparation: Sequencing generated READ_NUMBER 150000000 paired-end reads of READ_LENGTH 76 after amplification with PCR_PROBABILITY 0.05.Fragmentation used RNA with UR and FRAG_UR_ETA 350; reverse transcription was lossless with RT_MIN 500 and RT_MAX 5500.
Supplementary Figure 5: Convergence of relative abundance estimates
Supplementary Figure 5 shows that relative-abundance estimates converge rapidly during SQUAREM iterations in the Universal Human Reference tissue experiment, with the residual appearing converged before 30 iterations.
- Convergence of relative abundance estimates: Before 30 iterations, the SQUAREM residual appears to have converged in the Universal Human Reference tissue experiment.The figure plots the average difference between relative abundances estimated by successive EM-step applications against SQUAREM iterations, showing a rapid residual decline.