Question: Calculating TPM from featureCounts output
9
3.9 years ago by
nash.claire380
nash.claire380 wrote:

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??

rna-seq • 23k views
modified 10 months ago by janzrimec0 • written 3.9 years ago by nash.claire380

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.

13
3.9 years ago by
Kamil2.0k
Boston
Kamil2.0k wrote:

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.

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]

# 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)
return(tpm)
}
``````
1

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

Great Kamil! Thanks.

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

1

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.

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!

6
3.9 years ago by
Rob3.6k
United States
Rob3.6k wrote:

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

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.

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.