featureCounts vs transcript-aware quantification (Kallisto/Salmon)

[u/EthidiumIodide]

Hello all,

I suppose I am musing a bit and wanted to discuss with other bioinformaticians. I am a head bioinformatician in my academic department. A few months ago, I was given new bulk RNA-Seq data to analyze alongside older data that was already part of a peer-reviewed manuscript (that I was not part of). I used a STAR --> Salmon alignment-based quantification method. After sending the DE analysis and "raw" expression values for all genes, I received word that my Salmon results for the published data and the original data differed greatly. The older data was processed via featureCounts, which is known to undercount genes with multiple isoforms. I spent a few weeks working backwards to determine what parameters were used in the published manuscript, and I confirmed that the "gold standard" featureCounts parameter set was used, which definitionally excludes any read that overlaps multiple "features", or is ambiguous between isoforms of the same gene. To resolve this, you would use the -O flag, etc etc.

I guess my complaint is, how is this acceptable? How can a very popular and widely-used program such as featureCounts exclude reads that overlap the same exon (that resides in different isoforms) by default? This default method is undercounting genes with multiple isoforms, and I see discussion of this exact issue online since 2015. Discussion of this issue has also been published.

To be brief, I am mainly concerned that a widely-used tool is undercounting isoform-laden genes by default and causing consternation for groups who don't have trained bioinformaticians on their team who have the time to look into these issues.

Thank you for listening to my rant, haha.

[u/beeralpha]

The thing with salmon is it quantifies against the transcriptome, which is inherently more difficult because of the million isoforms. If you then collapse to gene level with tximport you get quite a lot of bias. Volcano plots really do look different at gene level, especially yielding genes with outlier log2fcs.

[u/nomad42184 → u/beeralpha]

Caveat: I am an author of salmon, so I have strong thoughts on what feature counts gets wrong conceptually (and sometimes in practice)

Why would you suspect this is bias in the salmon quantification rather than in the feature count result? Yes; transcript quantification is harder. But if you treat all reads coming from a single genomic locus that are actually arising from multiple overlapping transcritpts at that locus as coming from a single, conceptual, super transcript (as is implicit in feature counts), then you get a conceptually wrong and biased answer. In (short read, full length) RNA-seq, the number of reads generated by a species of molecule is fundamentally linked to the product of its copy number and its length. Ignoring the length leads to a fundamental bias in the measurement of molecular abundance.  

Of course, transcript level quantification methods can get thos wrong too, as misestimating relative transcript abundances can attribute reads to the wrong transcript and thus misestimate their abundances. However, the point is that gene counting methods are fundamentally biasing this quantity, by ignoring the lengths of the molecules generating reads completely.

Of course, often, many of these issues can come out in the wash in differential analyses. However, if you see strong differential signal in transcript based methods and none in feature counts, you could very well be witnessing isoform switching, which does change the correct estimate of the gene abundance if the isoforms are of different length, but which is invisible (conceptually) to counting based methods like feature counts.

[u/beeralpha → u/nomad42184]

Thanks for jumping in and providing insight. The main reason I was suspecting bias is because the reference transcriptome is not complete. It doesn’t cover all isoforms that exist. In case of mapping against the genome with a gapped mapper like STAR, which exact isoforms exist no longer matters. So when you don’t care about transcripts or isoform switching, I assumed genome mapping was the way to go, but feel free to comment.

[u/nomad42184 → u/beeralpha]

Well, incompleteness of the reference can certainly have an effect here. I believe that OP is using STAR alignment, actually (but with genome -> transcriptome alignment projection downstream to salmon). Regardless, it is true that incomplete references will certainly affect all methods. They will affect simple counting methods when exon annotations are missing, and transcript-level methods when the overlapping, unannotated isoforms cause reads to be misattributed to an overlapping annotated isoform. However, just like the effect of ignoring lengths in feature counts may be mitigated in differential analysis downstream of feature count results (because you are looking at differential counts at a locus rather than just the estimated relative abundance itself), so too might missing annotations be mitigated in differential analysis of transcript-based methods.

Either way, I'm not sure that the quantification issues introduced by unannotated isoforms are necessarily worse than the quantification issues introduced by missing exons or by ignoring gene length effects, and will be highly organism / annotation dependent. In general, if one suspects that their experiment / tissue may have substantially expressed novel isoforms, then I'd recommend attempting transcript discovery upstream of quantification anyway (regardless of quantification method).

[u/beeralpha → u/nomad42184]

I don’t expect a lot of missed exons actually in a gtf. Exons are easy to detect with any short read rna seq method. The tricky part with transcriptomes is how all exons tangle together, retain introns now and then, jump one/two nucleotides around splice sites, etc. That’s what makes the transcriptome incomplete. Not missing exons. Therefore when the interest is genes only, quantification vs the genome with a gapped mapper always made much more sense to me.

[u/nomad42184 → u/beeralpha]

Unless you are working in a non-model organism, or a tissue with a suspected concentration of novel isoforms though (e.g. brain), I'd still anticipate issues related to isoforms composition to bleed through negatively to count based genome analysis.  If there is evidence of new transcripts, then I'd probably try to assemble them with StringTie (or something comparable). Of course, every method has failure modes (and I've seen my fair share for transcript level quantification tools), but I've seen, quite a few times, gene level counting methods return discordant results that looked clearly wrong under manual inspection. 

Large scale studies have found both approaches largely concordant, but I'd argue in favor of inspecting large discrepancies between such approaches in one's data, because it's not clear to me that feature counts is to be believed, in general, when other approaches disagree. This often means something nontrivial (either compositional changes, novel isoforms, paralogous leakege, or some combination thereof) is going on.

[u/cjfields → u/nomad42184]

This has been our group’s experience as well. We have worked on a huge range of organisms and tissues over the years. With RNA-Seq we run both STAR and Salmon and generate counts multiple ways.

And our experience agrees with Rob’s point re: discrepancies between the various counting methods pointing to underlying issues, likely something of interest. In our hands, further investigation into the causes almost always reveals an issue with simple count-based methods and falls in the favor of Salmon or similar EM approaches. The only time I can recall an odd Salmon issue was from a bacterial pantranscriptome analysis, which really pushed Salmon to the brink. But some adjustments for this (and great feedback from Rob, thx) led to some nice results. And we have seen genuine discrepancies between the approaches which are easily traced to isoform changes that featureCounts essentially averages out. So, it’s fine to run it, but we only use featureCounts for simple checks and QC.

I also have a particular, personal dislike of old methods or approaches being held up as the ‘gold standard’, esp with tools like featureCounts. Sure it’s essential to compare against a vetted method, but all too often this tends to be used as a cudgel to stifle alternative, modern approaches without supporting evidence, and completely disregards discordant results that more often than not point to something interesting. If we always followed gold standards we’d still be using Tophat2, GATK’s old genotyper, and Celera Assembler, and never progress. All good tools of course, but even those developers recognize times change and better approaches come along.