# Calculating most abundant transcript from RNA-Seq data

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

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.

• I did not know about terminus, it looks really interesting. Could you explain why TMM/median ratio method are still inappropriate across samples? I thought they were designed for between sample normalization. – Tomas Bencomo May 14 at 16:11
• Don't get me wrong, I think both TMM and median-ratio are good enough for your approach, I just don't like to throw around the term normalization lightly. IMO normalization implies aligning dsitributions to allow for statistical testing, which those two don't fully do. If there is a huge disparity in the number and magnitude of expressed genes between conditions, or lots of zero inflation, those methods do not produce a clean normalization either. – Bastian Schiffthaler May 15 at 19:22

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.

• OP seeks to find the most highly-express transcript, not gene. The idea with the ranks though should be considered, and one could rank by the TPMs that salmon returns for each transcript. That avoids the use of between-sample normalization. So one would quantify with salmon, take the TPMs, rank transcripts per gene, and then use the transcript with the minimal aggregated ranks per between all samples. – ATpoint May 14 at 13:06