Dear List,

Our sequencing core recently switched from STAR to Kallisto. The genewise output was given as both TPM values and counts. Gordon Smyth has written that TPMs should not be used as input to EdgeR. In a recent paper (F1000 Research 2016,4:1521) Soneson, Love, and Robinson have presented evidence that the method they call simplesum_avextl for deriving reads-per-gene, from numbers-of-transcripts is as good or better for differential expression than alignment+featureCounts. In a comment on that paper, Lior Pachter advocates simply adding the reads assigned to each transcript to get the reads for a gene.

I would greatly appreciate Gordon's input, or that of someone from his group, as to whether there is a proper way to get counts from pseudoalignment-based numbers-or-transcripts as input to edgeR or limma-voom. Or alternatively, is it recommended that I generate the counts from the fastq files with an alignment program, such as Rsubread.

Thanks and best wishes, Rich Richard Friedman, Columbia University

“My father justifies mice. You need to make sure that you have enough mice for an experiment and that you do not have too many. He makes sure that no mouse dies in vain. You are not allowed to use chimps, so you have to use mice.”- Rose Friedman, age 22

Dear Rich,

I'll tell you what I use. I use:

1. Rsubread::align + featureCounts + limma or edgeR-QL for gene level expression analyses

2. Rsubread::subjunc + featureCounts + limma or edgeR for differential splicing analyses

3. 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.

I'm not in Gordon's group, but I'll still offer an answer:

Just follow the instructions outlined in the tximport vignette. They provide instructions on how to roll up transcript level counts into gene-level counts for your favorite DGE pipeline ... as long as your favorite is one of edgeR, DESeq2, or limma/voom.