4 months ago by
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
I'll tell you what I use. I use:
Rsubread::align + featureCounts + limma or edgeR-QL for gene level expression analyses
Rsubread::subjunc + featureCounts + limma or edgeR for differential splicing analyses
Salmon (with bootstrapping) + edgeR::catchSalmon for transcript level expression analysis.
The short answer to your question is that getting genewise counts from kallisto is ok as input to either limma and edgeR. Of course you should work from the counts and not the TPMs. Either the tximport genewise pipeline or else just adding up counts over transcript for each gene are both ok from this point of view.
In fact, getting total transcript counts from kallisto is in principle the same as running featureCounts on the same transcriptome annotation except that kallisto is restricted to perfect match alignments so featureCounts should manage to count more reads.
There's a much longer answer, but I don't really want to write at length here or now. I will say that kallisto makes too many assumptions and approximations for my taste and I personally wouldn't hesitate to use Rsubread instead. Just to name one issue of many, more than 10% of Ensembl transcriptome annotation is made up of duplicates, i.e., the exact same transcript sequence but with multiple entries under different names, and kallisto does no checking of this. Genewise inference is robust whereas transcript-level inference is inherently noisy, so using transcripts as an intermediary to get genewise counts seems a peculiar choice to me.
Responding to your comment about Soneson et al (2016), they were not actually able to show any improvement for tximport over featureCounts for any real dataset. Of course their main point is true. If a gene has a two transcripts of very different lengths, and the short transcript is expressed more highly than the long, and the two transcripts are both DE but in opposite directions, then it is possible for total genewise count and total gene expression (sum of transcript expressions) to show fold changes in opposite directions. My understanding is that the tximport refinement is intended to get the genewise counts to better align with sum of transcript expression. This scenario isn't very common in practice but if you run a simulation for which this scenario is the norm then of course you can show a difference. My view is that the tximport refinement isn't a complete solution, as it would be better to detect the discordant DE, and this refinement is in any case dwarfed by other issues with kallisto and imperfect annotation.