Calculating TPM from featureCounts output

Question

[nash.claire]

Hi all,

Have a simple question but just want to double check I'm not doing something stupid.

I have paired-end RNA-seq data for which I have used featureCounts to quantify raw counts. I now want to normalize using the TPM formula. I read this blog :-

http://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/

which simply says, divide read counts by gene length in kilobases to give reads per kilobase (RPK), sum all the RPK values and divide by a million for a per million scaling factor and then divide all RPK values by this scaling factor.

So taking my output from featureCounts which looks something like this : -

Geneid      Chr    Start          End            Strand      Length   sample.bam
NM_032291   chr1   66999639 etc   67000051 etc   + etc etc   10934    25

I use the values in the "Length" column as reads per kilobase or do I have to convert this to per kilobase first? It didn't say in the Subread manual much about this length value. And the value in my "sample.bam" column is definitely the read count value I need?

[janzrimec]

The paper to cite seems to be Bo Li et al. 2010: https://academic.oup.com/bioinformatics/article/26/4/493/243395 - see page 494 eq. 2, though it references a paper within, and the RNAseq review paper Conesa et al. 2016 references Pachter 2011: https://arxiv.org/pdf/1104.3889.pdf.

Answer 1

[Kamil]

Consider using the tximport R package:

https://bioconductor.org/packages/devel/bioc/vignettes/tximport/inst/doc/tximport.html

---

Below, I'm sharing a function I have used to convert counts to TPM.

See here for a script that gets the number of coding bases per gene.

See here for a script that gets the mean fragment length for each library.

[Shixiang]

What's the meaning of "meanFragmentLength", could you give an example or more description?

[Gautier Richard]

Fragment Length usually corresponds to the size of the DNA or RNA fragments that have been sequenced. Such information is usually obtained with a Bioanalyzer. It can also be infered from the distance between pairs when performing paired-end sequencing.

[Fernando Rossello]

Simple and clean, thanks for the function.

As a suggestion, shouldn't there be a parameter to keep feature names? As in:

counts_to_tpm <- function(counts, featureLength, meanFragmentLength, featureName)
# Ensure valid arguments.
stopifnot(length(featureLength) == nrow(counts))
stopifnot(length(meanFragmentLength) == ncol(counts))

# Compute effective lengths of features in each library.
 effLen <- do.call(cbind, lapply(1:ncol(counts), function(i) {
 featureLength - meanFragmentLength[i] + 1
}))

# Exclude genes with length less than the mean fragment length.
idx <- apply(effLen, 1, function(x) min(x) > 1)
counts <- counts[idx,]
effLen <- effLen[idx,]
featureLength <- featureLength[idx]
###  featureNameAdjusted <- featureName[idx]

# Process one column at a time.
tpm <- do.call(cbind, lapply(1:ncol(counts), function(i) {
rate = log(counts[,i]) - log(effLen[,i])
denom = log(sum(exp(rate)))
exp(rate - denom + log(1e6))
}))

# Copy the column names from the original matrix add feature names to rownames (or extra column?)
colnames(tpm) <- colnames(counts)
###rownames(tpm) <- featureNameAdjusted
return(tpm)
}

[Kamil]

Thanks for the suggestion! I modified the function for you.

[Fernando Rossello]

Great Kamil! Thanks.

[camillab.]

Hi so I am confused I want to calculate the TPM from features counts and you wrote a very nice script but do I have to calculate the mean fragment length necessary to get the TPM or can I just run instead:

tpm <- function(counts, lengths) {
  rate <- counts / lengths
  rate / sum(rate) * 1e6
}

[dganiewich]

Hi! From the link to the program for calculating fragment length I don't seem to find that exact metric. Instead, I found the read length and the insert length, do you use here insert as synonym of fragment? I understand that ideally only the insert is aligned and in that case they can be the same, but what happens if part of the adapter is also involved? I am a little confused about this, and why TPM accounts for the fragment in stead of the insert which is the actual DNA of interest... Can you give me a hand on this? Thanks!

Answer 2

[Rob]

First, let me suggest that you'd probably be better off using a tool for explicitly estimating relative abundance than processing the output of a tool like featureCounts (see e.g. this recent manuscript by Soneson et al.). To answer your question, that's a very round-about way of computing TPM, which seems to introduce some arbitrary scaling factors for no real reason (e.g. why compute RPK rather than just reads per base?). So, anyway, you can compute TPM in the following manner. Let N_i be the number of reads coming from gene i and L_i be the length of gene i. Then TPM_i is given by the following formula:

TPM_i = 10^6 * (N_i / L_i) / [Σ_j (N_j / L_j)]

That is, you first compute the rate of fragments per base for each gene (given by N_i / L_i for gene i), then you divide by the sum of this quantity for all genes --- this normalizes the measurements so they sum to 1. Now, you have a measurement with a non-negative value for each gene that sums to 1 (is a partition of unity); this is called the "transcript fraction" (or, in your case, "gene fraction"). You simply multiply this by 1 million to convert it to "transcripts per million" (TPM).

[CC]

Thank you for this explanation of how to calculate TPM. I was previously using the method described on the website linked in the question above - it gave me quite a different answer. I was wondering - is there a publication which lists the method to calculate TPM in the way you have in your answer? I've had a look but keep finding different methods to calculate it. Thanks for letting me know.

[enxxx23]

It looks like CPM is computed here.

TPM for a gene should be computed using the TPMs computed from transcripts.

It looks like a big confusion is going on regarding CPM, TPM for genes, and TPM for transcripts.

Answer 3

[gernophil]

Is there any way, how I can calculate TPM from count data from a publicly available dataset, if I don't have the meanFragmentLength? I can't find it anywhere.