From RNA-seq raw counts to normalized data

Question

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!

Answer 1

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?