1
$\begingroup$

I have several RNA-seq datasets. Some of them provide RNA-seq raw counts, some provide FPKM, RPKM and some have transcripts per million (TPM) data. Non of them provide fastq files, all data is processed already. At the end I want all datasets to be normalized to TPM.

I'm using this code in order to normalize raw counts to TPM: (using R)

rpkm <- apply(X = subset(dataset),
                MARGIN = 2,
                FUN = function(x) {
                  10^9 * x / genelength / sum(as.numeric(x))
                })
  
TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

When RPKM is provided, and no raw counts is available, I use the second line in the same code:

TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

And when FPKM is provided, I use this formula to transform the data to TPM:

TPM = FPKM*X 

where X = 1e6/[sum of all FPKM of a sample]

NOTE: genelength is obtained using the biomart package in R, to get the transcript length directly out of ensemble.

However, because I know the steps of TPM normalization in theory, one should firstly normalize to gene length, and then to gene depth. I'm not sure this is what the code does, feeling very skiptical about it even thought this code was given to me by a Phd student.

Can you guys please help me with this question? is my code right or should I alter something?

Thank you!

$\endgroup$
1
  • $\begingroup$ Regarding the use of TPM, it's a good idea to be aware of the issues associated with it, even if it's your only available option. There are a few different TPM questions that have been asked previously, this answer has a reasonable summary. $\endgroup$
    – gringer
    Oct 15, 2022 at 9:11

1 Answer 1

2
$\begingroup$

However, because I know the steps of TPM normalization in theory, one should firstly normalize to gene length, and then to gene depth.

No, you should normalize to transcript length. Which requires you to know the relative abundance of all the transcripts so you can correct for all their lengths.

I do not think it's possible to un-normalize RPKM or TPM without knowing all the lengths involved in doing the calculations.

Let me give an example: here is RSEM output for just one gene. Note that the software thinks there are three transcripts present for this one gene:

transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct

ENST00000373020 ENSG00000000003 3768    3718.54 179.95  13.88   7.94    48.61

ENST00000494424 ENSG00000000003 820 770.54  0.00    0.00    0.00    0.00

ENST00000496771 ENSG00000000003 1025    975.54  11.94   3.51    2.01    12.29

ENST00000612152 ENSG00000000003 3796    3746.54 0.00    0.00    0.00    0.00

ENST00000614008 ENSG00000000003 900 850.54  33.11   11.17   6.39    39.10

ENST00000373031 ENSG00000000005 1205    1155.54 0.00    0.00    0.00    0.00

ENST00000485971 ENSG00000000005 542 492.54  0.00    0.00    0.00    0.00

So, what would you say the effective length should be for the gene level calculations, and would you get that right without having this transcript-level breakdown?

$\endgroup$
5
  • $\begingroup$ I do have the transcript lengths, as I mentioned above I got them from ensemble. $\endgroup$
    – Kev
    Oct 18, 2022 at 8:08
  • $\begingroup$ Do you have transcript level counts? $\endgroup$
    – swbarnes2
    Oct 19, 2022 at 4:56
  • $\begingroup$ No.. I assume this is something you get when you have BAM files? Anyway, I don't have that. I was able to get the transcript length which is necessary for the normalization. Do you think the code here provided by Michael Love is better? support.bioconductor.org/p/91218/#91256 $\endgroup$
    – Kev
    Oct 19, 2022 at 7:17
  • $\begingroup$ No, bam files do not magically give you transcript level counts. Software like Kallisto and RSEM and Salmon do, and they do all the calculations for you. Transcript length is not enough. Genes have multiple transcripts of different sizes, and you can have a mix of isoforms in your samples. The software I mentioned will magically handle all this for you. I don't know that it's possible for you to do the magic in reverse. I'm sure you can wave your hands and do some math and get numbers. I don't know that those numbers will be accurate. $\endgroup$
    – swbarnes2
    Oct 19, 2022 at 16:57
  • $\begingroup$ I see. According to Michael Love, as you said, there is no way to calculate accurate TPM values out of count matrix. As I mentioned earlier I have no access to BAM files unforunatlely. But in this link I think I got my answer, they discuss this matter, and Michael Love provides here a simple code to normalize TPM out of count matrix. Thank you for your help and explenations. This is the link support.bioconductor.org/p/91218/#91256 $\endgroup$
    – Kev
    Oct 20, 2022 at 10:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.