Question: Calculating FPKM after htseq-count
1
gravatar for Fill
2.6 years ago by
Fill70
Fill70 wrote:

I just want to make it clear.

I need to calculate FPKM. I use this formula: Normalized = [(raw_read_count)(10^9)] / [(gene_length)(XXXX)],

XXXX = the count of all reads that are aligned to protein-coding genes in that alignment.

How should I calculate XXXX? Is it just sum of all raw_read_counts after htseq-count (e.g. in R it will be XXXX <- sum(collumn_with_raw_red_counts)?

Thanks!

rna-seq fpkm htseq • 8.2k views
ADD COMMENTlink modified 5 months ago by Ahmed Alhendi30 • written 2.6 years ago by Fill70
1

Try this solution,

PERL solution:

https://github.com/santhilalsubhash/rpkm_rnaseq_count

R solution:

A: How to normalise read count per gene

ADD REPLYlink modified 2.6 years ago • written 2.6 years ago by EagleEye6.4k

Thanks, but the question is about how to get Total count of all reads. Is it just sum of all counts? I ask it because I am trying to duplicate the results on GDC data portal and I can't do it because their Total count of all reads smaller than mine (which I calculate by sum())

ADD REPLYlink written 2.6 years ago by Fill70

Yes. Sum of all counts from individual samples (library size).

ADD REPLYlink written 2.6 years ago by EagleEye6.4k

It's "fragments per kilobase transcript per million reads", so you should divide by million, not multiply.

That said, are you sure FPKMs is something you need? I don't know your application, but for many purposes, there are better normalisation methods.

ADD REPLYlink written 2.6 years ago by WouterDeCoster40k
1

formula This is formula with effective length, you can see why I multiply by 10^9, not devide.

I am trying to duplicate the results on GDC data portal.

ADD REPLYlink modified 2.6 years ago • written 2.6 years ago by Fill70

Ugh, you're totally right. Shame on me!

ADD REPLYlink written 2.6 years ago by WouterDeCoster40k

Can you clarify duplicate part? Unless you are using the exact versions of software/genome build that may not be realistically possible.

ADD REPLYlink written 2.6 years ago by genomax71k

I am using exact versions software and genome like TCGA. And everything is good (I mean results from TCGA and my results after mRNA pipeline are identic. Except for N value for calculating FPKM (the count of all reads that are aligned to protein-coding genes in that alignment)

ADD REPLYlink written 2.6 years ago by Fill70

Is it possible that for paired end reads (fragments) they divided total reads by 2 ? (I know it is silly to think like that but can you match your numbers divided by two with TCGA total library sizes)

OR

While calculating FPKM (raw_reads/2) per gene.

ADD REPLYlink modified 2.6 years ago • written 2.6 years ago by EagleEye6.4k

Thanks for your guess, but it doesn't work that way. My N values differs by 0.03x - 0.08x (x = TCGA N values).

ADD REPLYlink modified 2.6 years ago • written 2.6 years ago by Fill70

Aligners may produce non-deterministic output (unless they are able to accept a seed). Perhaps that is what is causing this difference.

ADD REPLYlink written 2.6 years ago by genomax71k
6
gravatar for Fill
2.6 years ago by
Fill70
Fill70 wrote:

I've got answer from GDC portal:

  1. Download GTF files used in HTSeq analyses: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files (GDC.h38 GENCODE v22 GTF)
  2. Extract only protein-coding gene IDs: less gencode.v22.annotation.gtf | grep "\tgene\t" | grep protein_coding | cut -f9 | cut -f2 -d '"' > EnsembleIDsPCG.txt
  3. Use resulting list to extract only protein-coding values from counts file: less CountFile.txt | grep -Ff ProteinCodingGeneList.txt > CountOnlyProt.txt
  4. Sum the values of "CountOnlyProt.txt" and that will give you your denominator value.

My problem was that I counted reads for all genes, but should only for protein-coding.

P.S. thanks to GDC support team!

ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by Fill70

Hi, I have stumbled upon your post trying to find out why I am not able to obtain the FPKM values provided by the GDC using the same raw count data. After following these steps it does not get any better. Just to be sure, you took the counts from the HTSeq files and the gene lengths from the GDC.h38 GENCODE v22 GTF, right? As for the N in the denominator, as explained above, it should be the sum of all protein coding genes...

ADD REPLYlink written 23 months ago by CuriousGuy30

You're right. Sum of all reads which align protein-coding genes. Gene lengths and counts are for all genes.

ADD REPLYlink written 21 months ago by Fill70

Does it make sense to also include genes that code for non-coding regulatory rna? To calculate N by summing counts in all exonic features?

ADD REPLYlink written 4 months ago by gatollefson0
0
gravatar for Ahmed Alhendi
5 months ago by
University of Southampton, UK
Ahmed Alhendi30 wrote:

Try countToFPKM package. This package provides an easy to use function to convert the read count matrix into FPKM values normalised by library size and feature effective length. Implements the following equation:

enter image description here.

The fpkm() function requires three inputs to return FPKM as numeric matrix normalized by library size and feature length:

  • counts A numeric matrix of raw feature counts.
  • featureLength A numeric vector with feature lengths that can be obtained using biomaRt.
  • meanFragmentLength A numeric vector with mean fragment lengths, which can be calculate with
    Picard using CollectInsertSizeMetrics.

Also see https://github.com/AAlhendi1707/countToFPKM

ADD COMMENTlink modified 5 months ago • written 5 months ago by Ahmed Alhendi30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1993 users visited in the last hour