I tried using edgeR::rpkm() on the feature counts file like this post (unfortunately the link you supply there to your RNA case study isn't working for me)
#read in input feature counts file with ENSGID, counts, and gene length columns
file <- read.delim("D0-Fibroblast_rep1 GSM4553868_19-03967_CCTAAGAC.txt", header = TRUE, sep = "\t", row.names = 1)
#create DGEList object containing the counts and gene IDs as row names
counts <- DGEList(counts=file$X19.03967_CCTAAGAC, genes=rownames(file))
#Calculate scaling factors to convert raw library sizes into effective library sizes
counts <- calcNormFactors(counts)
#calculate rpkm from counts in DGEList object
rpkm <- rpkm(counts, gene.length=file$geneLength)
#add geneIDs as rownames to rpkm object
rownames(rpkm) <- rownames(file)
Is this the correct way to use the rpkm function? And would the output essentially be FPKMs?
For context, downstream I'm trying to use this RNA data for Paired Chromatin Expression Analysis which states it needs FPKM/TPM values. Thank you for your help.
• 0 views
•
link
written
by
cthangav •
I'm not sure how those lengths are determined. Ensembl doesn't think any transcripts of that first gene are 4kb long.
It's gene length, not transcript length. It is the total number of bases covered by any annotated exon for that gene. It merges all the annotated transcripts.
But would you want to use that value for TPM? Or anything?
It works well for RPKM or FPKM. It is the appropriate length measure here because it matches how the genewise read counts are computed.
What must be clear to everyone using RPKM is that simple summing annotated exons is a comparably naive way of getting the gene length as it does not respect in which proportions isoforms are being expressed. Hence, if you want to make statements on expression levels with RPKM this is not accurate. The concern is not theoretical, for example the transcription factor CEBPa (a known driver in hematopoietic malignency, and the isoforms having known and different roles in it) expresses two isoforms with ~25% difference in tx length. It then makes quite a difference whether conditionA expresses 95% isoform-1 and 5% isoform-B but conditionB expresses a 10/90% ratio. featureCounts (please correct if wrong) in my understanding will output the same gene length for all conditions. If this unequal ratio is nested by experimental conditions then RPKMs on naive exon sum is not accurate and one needs dedicated approaches like salmon or rsem to estimate effective length based on transcript expression composition. I would therefore treat featureCounts length output with care if @OP you really want to use that. Whether this affects the big picture is questionable and would need to be checked, it definitely can affect individual genes.
Yes, I am aware of the approximations made by FPKM. (I've been aware of them for 15 years, although people do love to explain it to me again.)
IMO FPKM will probably work fine for OP, just as it did for the original authors of the PECA method that OP is attempting to reproduce. The PECA method is from Wing Wong's group, which is a very strong statistical bioinformatics Lab at Stanford.
This is not the time for a long discussion, but I do wish that people who answer questions on Biostars would treat TPM with a similar level of scepticism. In my experience, FPKM seldom gives a misleading idea of expression changes for more than 1-2 genes per dataset. On the other hand, transcript TPM methods like RSEM or Salmon give TPMs that are wildly wrong for tens of thousands of transcripts for almost every dataset. That is not a criticism of those methods, just an observation on what can be reliably measured from available data.
Dear Professor Smyth,
Can you show me referenced benchmarks that FPKM is misleading for "1-2 genes per dataset" while TPMs "are wildly wrong for tens of thousands of transcripts for almost every dataset".
There are many papers (w/ benchmarks) arguing the limitations of FPKMs and the union-exon approach to gene length (and why TPM is preferable) -- Zhao et al., 2015, the tximport paper, Trapnell et. al 2012, and the blog post (w/ included CSHL talk) from the person who first formulated FPKM which, as you are aware, is a "statistical bioinformatics lab" at Berkeley (now Caltech).
Hello Ram and ATPoint, thank you for helping me repeatedly with my questions. I am now finalizing the paper, and if there is a way to acknowledge you, I would be happy to do so. Thanks again!
Best way to get TPMs is to use RSEM instead of featureCounts.
The notion of gene length is kinda iffy. Does it mean the length of the entire gene (introns+exons)? Does it mean the lengths of exons? What if one transcript of a gene is expressed very highly while a shorter one of the same gene is not expressed at all -- do we account for that?