Calculating most abundant transcript from RNA-Seq data

Question

vcf2maf uses VEP to annotate variants, and I believe selects the default Ensembl transcript to use for annotation. Sometimes the transcript that VEP selects is not the transcript I'm interested in, usually because the selected transcript is not the most highly expressed transcript in my tissue of interest (skin). vcf2maf allows you to provide a transcript override list so that VEP annotates the variant using the specified transcripts instead.

I have several skin samples sequenced with bulk RNA-Seq. I want to estimate the average abundance for each transcript across all samples and then use these abundances to rank transcripts from most to least abundant. Then I will use the most abundant transcript as the default VEP transcript. I plan to use salmon or kallisto to quantify transcript abundance. Should I use TPM or normalized counts to calculate average expression?

My initial thought is to use normalized counts (generated by DESeq2 from raw counts). Are there any problems with this approach? GTEx displays transcript abundance with average TPM, but I thought TPM was inappropriate to use across samples because it doesn't account for between sample differences.

Update: I forgot to mention I also tried using TPM ranks like @ATpoint describes. I haven't fully compared how this compares to transcripts identified by normalized counts, but the initial genes I checked showed good concordance between methods

Answer 1

Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here:

1. salmon (as you mentioned) to quanitfy. Run it with --numGibbsSamples 50 (or higher if your computer has the hardware. Also run with -d for the next tool
2. terminus is a tool that also originates from the COMBINE-lab, It's relatively experimental, but it applies well to your use case. It will take the salmon output and assign counts to transcript groups. If a set of transcripts can be clearly quantified above a certain confidence, that transcript will stay unique. If it is confounded with other transcripts, those will be placed in a transcript group and their counts will be summarised. More info here

You are right in that TPM, FPKM etc. are not appropriate measures to compare across samples. TMM, or DESeq2's median ratio method (what you get with counts(dds, normalized=TRUE) are a step in the right direction, but are still inappropriate across samples. I'd suggest this often posted blog post in the topic.

There are no perfect solutions here since RNA-Seq is inherently a relative quanitifcation method so this will always be an issue. IMO, use DESeq2::varianceStabilizingTransformation(), DESeq2::rlog() or limma::voom() after terminus and rank your transcripts that way.

Just keep the inherent transcript incertainties in mind.

Answer 2

If I read you correctly, you want to find the gene with the highest expression?

With RNA-seq data from multiple samples, I'd find this by ranking each gene by TPM/FPKM and taking the mean rank across samples. This method will tend to mitigate any large differences in a gene's expression between samples.

As you are comparing (RNA-seq) expression across genes you must take into account transcript length.