如何计算cuffdiff中的FPKM值 65

FPKM, 是expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced缩写。直译过来就是每百万测序碱基中每千个转录子测序碱基中所包含的测序片断数。与RPKM不同的是,RPKM是直接使用的reads数,而对于FPKM,如果是pair-end的话有可能有些mapped至基因组的是一对,那就算一个片断,如果map至基因组的是只有一侧的read,也算一个片断。

FPKM的公式就可以从字面上写成:10^6*10^3*\frac{C}{NL},其中C是map至该基因的外显子上的片断数,N是所有map至基因组的测序reads的碱基数,L就是该基因外显子碱基全长。

在其文献(Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation)Supplementary Text and Figures中是这样描述的:
FPKM=10^6*10^3*\beta _g*\sum_{t \in x} \frac{\bar{\gamma _t}}{\grave{l_t}}
其中beta和gamma都是likelihood function中的参数。lt被定义为\grave{l_t}=\sum_{i=1}^{l(t)}F(i)(l(t)-i+1)因为,假设transcript t的长度为l(t),那么随机分布时在某一位置t出现长度为k的片段的概率为\frac{1}{l(t)-k}

在cuffdiff中,它会将同一组中所有的样品试为同一来源样品处理,这就是为什么同一组内如果有三个样品的话,最终得到FPKM值并不是三个样品单独运算得到的FPKM值的平均值。

65 thoughts on “如何计算cuffdiff中的FPKM值

  1. Reply Bo 4月 13,2013 10:47 上午

    “在cuffdiff中,它会将同一组中所有的样品试为同一来源样品处理,这就是为什么同一组内如果有三个样品的话,最终得到FPKM值并不是三个样品单独运算得到的FPKM值的平均值。”
    您好,我曾经试着分开运算每个样品,确实发现总的结果不是平均值。可是如果将所有样品作为一个sample来处理,那后面的p值是怎么得到的呢?我的数据里有些基因在一个sample里表达量是0,但在同一组另一个sample里FPKM却有几百。不知道这个是否正常?

    • Reply admin 4月 14,2013 8:58 下午

      你问到有些基因在一个sample中的表达量与其它同组的sample的表达量差别很大,是否正常。这完全有可能是正常的。我时常遇到这样的数据,但最后得到的结果验证都还不错。这就是生物样品的复杂性,即使平行实验,也可能差异很大。

  2. Reply Bo 4月 13,2013 10:55 上午

    不好意思,我还有一个问题。就是在运行玩cuffdiff以后发现对应FPKM最高值的大部分基因居然都没有gene id。我的预期是高表达量的区域应该对应着house keeping gene,那么怎么会还没有id呢?另外就是没有gene id的转录区域是不是就是还未发现的新的转录区呢?谢谢!

    • Reply admin 4月 14,2013 8:55 下午

      这个问题比较复杂。一个是检查你的注释文件是否正确,是否在cuffdiff时正确使用。第二就是查看那些没有gene id的转录区域是否真的没有已知注释。但是没有gene id的转录区并不意味着一定就是新的未知转录区。

      • Reply Bo 4月 15,2013 2:51 下午

        谢谢你的回答!关于注释文件是否正确,确实,在我run cuffmerge和cuffdiff的时候都出现了“warning: Couldn’t find fasta record for ‘chr13_random’”…还有好几个chromosome(但不是所有)也有一样的问题。我试着从不同的地方下reference,但结果都是这样。我在网上看到有人和我遇到了相似的问题,但都没有好办法解决。不知道这个问题会不会有影响。另外我发现之前我在问题里写错了,不是没有gene id,而是没有gene name

  3. Reply Bo 4月 16,2013 9:35 上午

    哦。谢谢啦!那么关于我的另一个问题,就是既然cuffdiff里的value具体是怎么算的呢?当我对每个sample单独运行cuffdiff时,发现个别基因在整个运行时得出的value比每一个sample单独运行cuffdiff时的value都要小,说明总的value确实不是平均值。 还有那个p值是基于每个sample的FPKM算出来的吗?

    • Reply admin 4月 17,2013 10:08 下午

      在其说明中这样写的:To calculate a p-value of observing the real log-fold-change under this null model, we simply sort all the samples and count how many of them are more extreme than the log fold change we actually saw in the real data. This number divided by the total number of draws is our estimate for the p-value.

  4. Reply Bo 4月 18,2013 10:14 上午

    好的。谢谢你的回答。

  5. Pingback: 如何统计BAM文件中的reads数 ← 糗世界

  6. Reply zhe 1月 18,2014 4:34 下午

    如何能知道cuffdiff使用了多少tophat结果中的reads呢?也就是说,比对的结果中有多少reads对于计算差异表达是useful的呢?也就是说,如何count有多少reads能比对到cuffdiff所使用的GTF上?

  7. Reply zhe 1月 22,2014 10:43 下午

    谢谢admin!

    但是我发现了一个现象,对于A、B两个样品,针对特定一个GTF,有时候,A比对上的fragment比B的多(多1%),但分别对A和B的FPKM求和后,居然发现A的FPKM之和小于B(少了近20%),不知您怎么看这种现象?

    • Reply admin 1月 23,2014 9:01 上午

      有很多种原因可能导致这一结果。1.非unique mapped reads。2.single ends数多。3.cufflinks发现了新的transcripts。我想到的暂时就这么多了。

  8. Reply gm 3月 11,2014 5:42 上午

    您好,有个问题想请教下。
    就是cuffdiff的作用到底是干嘛的,它不是进行sam和sam之间的比对吗,为什么在输入文件中还要提供一个.gtf文件?这个.gtf文件起什么作用呢?
    我是个小白,还希望您不吝赐教。谢谢。

    • Reply admin 3月 14,2014 3:05 下午

      对于map好的测序结果,一般只包括reads在基因组中的位置,并没有注释到每个基因有多少个reads。在比较时,输入gtf文件让cuffdiff可以依据一定的模型来计算对应到每个具体的基因上会有多少个reads。

  9. Reply orange 3月 12,2014 1:59 上午

    您好,我想问下,可以直接从merged.gtf文件中利用什么工具计算出FPKM的值吗?或者是,前面使用cufflink后的结果文件transcript.gtf中是有FPKM的,但这要怎么一一对应到merged.gtf文件中去呢?其实我现在思路也有点乱,不知道该怎么问比较合适,或许您能理解我的意思,麻烦了!谢谢~

    • Reply admin 3月 14,2014 3:10 下午

      你可以使用htseq-count来获得一个计数表。如果需要FPKM值的话,你可以在cuffdiff的结果中找到,有一些track文件,其中有每个样品的FPKM值。

  10. Reply qin 3月 16,2014 8:57 下午

    你好,我想问下,用cuffdiff分析差异表达基因时,如果有4个样本,两个病例,两个对照,分析时的merged.dtf是4个样本合在一起的转录本么?要比较病例和对照之间的差异表达,cuffdiff语句后面跟的是4个样本的.bam文件么?

    • Reply admin 3月 17,2014 8:03 上午

      看你的需要。如果你需要发现新的转录本的话,是需要4个样本合在一起的转录本。如果是其它情况,有可能只需要已知的转录本。在比较时,应该把所有涉及到的bam文件都传进去。

  11. Reply qin 3月 17,2014 11:46 下午

    谢谢,我还想问下,4个样本中2个对照,两个病例,要找出他们的差异表达基因,但cuffdiff运行得出来的结果gene_exp.diff文件中是分别作了4个样本间两两的差异表达,怎样可以得到病例组和对照组的差异表达基因啊

  12. Reply qin 3月 18,2014 7:18 下午

    不好意思,我还想问下,这样做出来的差异表达基因以Q值<0.05和log2(fold change)绝对值大于2删选,结果用做通路富集,可是没有富集的通路是怎么回事啊

  13. Reply qin 3月 25,2014 1:47 上午

    谢谢,如果要找出所有的差异表达基因,包括已知的基因和异形体,新的异形体,是不是要分两种步骤,一是利用cuffdiff有参考基因组的找出已知的异形体和已知的基因,另外一种是不是要利用cuffcompare过滤掉已知的,再进行cuffdiff分析?第二种的结果里面全是新的吧?

  14. Reply lexb 4月 11,2014 3:45 上午

    太感谢分享了,我想问个问题。cufflinks在计算FPKM的时候,如果read1比对上chr1,read2比对上chr9,这样这么算呢?
    看了比对结果统计,感觉Reads mapped in proper pairs不是很高,这么多这样的情况cufflinks是怎么算的啊?

    • Reply admin 4月 11,2014 2:29 下午

      paired end在map的时候就会要求输入两个reads之间可能的距离,所以并不会出现你所谓的read1在一个chromosome,而reads2在另一个chromosome。

  15. Reply lexb 4月 21,2014 8:47 下午

    您好,请问一下链特异性建库怎么找反义转录本。cufflinks可以做吗?

  16. Reply ziliu 5月 7,2014 7:37 上午

    您好!我是用tophat-cufflinks-cuffdiff 这个流程来做差异表达,但是因为样本数太多(cancer VS normal 70对),cuffdiff 运行得非常慢。现在有几个问题想请教下:
    1.cufflinks 的结果中除了cuffdiff需要的gtf,还是有每个gene 和isoform 的FPKM,为什么cuffdiff 还要从bam 文件出发,计算差异表达基因?
    2.FPKM是不是只相当于做了样本内的normalization,而算差异表达基因是需要样本间normalization的
    3.有没有什么方法可以利用每个样本cufflink FPKM的值 自行计算差异表达,达到cuffdiff的效果?

    • Reply admin 5月 7,2014 8:08 上午

      对于样品特别多,请使用cluster或者云计算。
      1. cufflinks在设计之初需要考虑两种情况,一种是发现新的转录本,一种是不需要发现新的转录本。如果直接是第二种,可能就不会经过cufflinks以及cuffmerge,所以在参数设置时,不可能假设你已经有了这些输入。第二,cufflinks给出的结果没有经过cuffmerge这一步,不具有比较的意义。
      2. 正确。但是,FPKM本身是会依据cuffmerge的结果而变化的。有时候因为比对的样品不同,对于相同的样品FPKM值也会不同(当然差别不是很大),你做得多了会慢慢体会。所以,你可以这么理解,但又不全然如此。
      3. 你可以查看htseq-count + R/bioconductor:::edgeR(or DESeq)。但是没有cluster的帮助,你计算量很大的情况下,一定非常慢。即使有cluster,你也要把心态放宽,1,2个星期是正常的,我时常要等1个月才能看到结果。

  17. Reply ziliu 5月 7,2014 8:12 下午

    非常感谢!不过我还有个问题不明白:
    假如我不考虑新的转录本。可以用cufflinks 的FPKM结果,自行做样本间normalization 和差异表达检验。这个与cuffdiff出来的结果差别吗?
    因为HTseq-DESeq这个流程只能给出差异表达基因,而cuffdiff还能给出差异表达的isoform,所以我特别想利用cufflinks isoform.tracks 这个结果,但是没有计算资源支持我用cuffdiff。

  18. Reply ziliu 5月 7,2014 8:13 下午

    sorry 这是个问句 假如我不考虑新的转录本。可以用cufflinks 的FPKM结果,自行做样本间normalization 和差异表达检验吗?

    • Reply admin 5月 8,2014 8:25 上午

      你可以考虑使用DEXseq。这个可以得到差异表达的isoform,不过可能会比较复杂。而且在性能上也不一定会比cuffdiff强。还有,如果使用DESeq,还是使用raw reads count会好些。这些信息可以在track file里得到。

  19. Reply Yamol 5月 15,2014 2:45 上午

    博主您好,我现在有一个treat,里面重复得到2个结果;同样control也有2个结果。我先用tophat跑完之后,用cufflinks分别跑了结果,得到4个文件。现在我该怎么做cuffdiff呢?

  20. Reply jiyang 6月 10,2014 6:54 下午

    我想问两个问题:(1)如何用cufflink软件确定一个新基因。具体的code或是从哪些结果中去分析。
    (2)如何单独计算一个没有gene_name的基因的FPKM。

  21. Reply Li 7月 9,2014 12:55 下午

    您好,请问
    1、用cuffdiff计算出来的gene_exp.diff中每一个基因的value值是不是就是FPKM?
    2、“在cuffdiff中,它会将同一组中所有的样品试为同一来源样品处理,这就是为什么同一组内如果有三个样品的话,最终得到FPKM值并不是三个样品单独运算得到的FPKM值的平均值。”
    同一来源样品处理,是不是就是说把3个重复的合并成一个样品,合并之后再计算FPKM,这样的话p-value又该怎么算?还有一个q-value是什么呢?
    3、虽然我看了您在回答中关于p-value的一段引用,对于它的计算还是觉得挺模糊,请指教。

    • Reply admin 7月 9,2014 3:52 下午

      1,是FPKM。
      2,你可以查看gene.fpkm_tracking文件。
      3, 这个问题比较复杂。你可以查看http://en.wikipedia.org/wiki/Jensen%E2%80%93Shannon_divergence和http://en.wikipedia.org/wiki/Delta_method来了解更多。其实我也没在太搞清楚。有机会再交流一下。
      4, significant的判断是依据FDR的值。The default is 0.05.

  22. Reply Li 7月 9,2014 12:57 下午

    另外对于结果中的significant的yes/no又是怎样判定的?

  23. Reply liu 7月 17,2014 2:09 上午

    您好,我也是用的tophat-cufflinks-cuffdiff这个流程,请问一下,cuffdiff结果里一些转录本的FPKM值都为0可能是什么原因?

  24. Reply mihuzx 7月 30,2014 8:20 上午

    博主:
    您好,
    感谢您的这么好的技术博客,学习了很多。
    现在有个问题:我用cufflinks 来计算FPKM值,但结果和Raw Count(由htseq-count计算) 出入很大,比如:四个样本的RawCount分别是0,0,5744,4606 而RPKM却是0,0,0,402.494。这出入也太大了啊,请教我该怎么办。重复了几次都是这样,但也没见报错。
    还有博主还有没有其他计算RPKM和FPKM的工具,给推荐一下,非常感谢。

    • Reply admin 7月 30,2014 12:05 下午

      其它的工具我也没有了,主流就是这两个。RPKM值和rawCount之间的转换的确挺难的。如果你给的值如果不是0多输入了一个的话,我觉得是非常正常的结果,没有出入很大啊。你需要正确理解RPKM,它和raw count不完全一样,它是标准化之后的值。

  25. Reply weikf 8月 3,2014 11:59 下午

    我最近用cuffdiff,碰到个问题。

    输入的gtf文件中某个geneid的区间是:
    chr16 BED exon 93792396 93793627
    chr16 BED exon 93794815 93795133
    但是cuffdiff的结果:

    区间是:
    NONMMUG017219 chr16:93696959-93810833

    请问,这是怎么回事呢?

  26. Reply xuanhe 8月 15,2014 2:23 上午

    博主您好:
    我想问一下cuffdiff的结果中fold change为多少的时候可以确定基因在一个样本中是特异的?还有cuffdiff结果中的P value是怎么算的呢?泊松分布?另外,如果在多个样本(比如6个)之间筛选特异的基因该怎么做呢?是将两两比较的结果合并到一起再筛选?这里面应该有一个normalization的问题吧?十分感谢!

    • Reply admin 8月 15,2014 9:37 下午

      p值的计算比你想象的要复杂,不是简单的泊松分布就可以的。normalization的过程在最初的FPKM值的时候就已经做过了。具体请阅读原始文献。

  27. Reply qianmacao 12月 6,2014 8:53 下午

    您好,我想请问一下我按照nature protocol上的那篇文章的步骤用tophat-cufflinks-cuffdiff这个流程来找差异表达的基因,但是最后得到的结果gene_exp.diff中significant没有为yes的,可能有什么原因呢?
    谢谢博主!

    • Reply admin 12月 8,2014 8:57 上午

      原因有很多种,最可能出现的问题在于平行样品之间差异过大,或者比对样品之间差异过小。其它原因,比如reads数不足,或者说基因组覆盖度不够,或者你的脚本中有小错误等等。

  28. Reply qianmacao 12月 18,2014 4:11 上午

    您说的那些原因应该都不是,因为数据还有实验步骤我都是按照《Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks》这篇文章做的,与原文数据结果比较后发现,我得到的gene_exp.diff文件中log2(fold_change)中有很多inf和-inf的值,但作者的数据中都是很大的数值,不知道为什么?
    还想请教您一个问题:Cufflinks是具体通过什么来判定一个基因为差异性表达的基因(significant==yes)的呢?谢谢!

    • Reply admin 12月 19,2014 2:41 下午

      它先通过 Poisson model来计算出一个pvalue,然后用这个pvalue较正一下得到qvalue,也就是FDR,然后使用fold change 和qvalue来判断是否significant.

  29. Reply shigongyao 12月 31,2014 3:06 上午

    博主新年好!
    我有个问题想请教您。我有棉花两个亚种A2和D5的基因组信息,知道它们之间转录本的对应关系,现在有它们分别10天的纤维转录组数据。怎么比较它们两者间的差异。这个好像不同于一个参考基因组多个样本的实验,而是两个参考基因组对应两个样本。怎么操作啊?cuffmerge和cuffdiff怎么做?
    谢谢!

  30. Reply zhihuawu 2月 15,2015 2:02 上午

    博主,你好!
    cuffdiff结果如何做样本相关性重复分析?还有cuffdiff结果不会显示每一个重复的表达倍数变化,有没有可以计算每个重复样本的表达变化?谢谢

  31. Reply Zhonghua_li 6月 11,2016 10:56 下午

    楼主您好,我的问题是两组比较中同一样品中的同一基因的表达量不同,但是相差不大,这是什么原因呢,是cuffdiff在比较之前有预处理吗?

  32. Reply 若凡 6月 21,2016 8:50 下午

    楼主您好,非常感谢您的博客,学到了很多东西。
    关于FPKM,我一直都有个疑问:FPKM计算得到的是每个转录本(isoform)的表达量,一个基因可以有多个转录本,这个基因的每个转录本(isoform)分别有个FPKM,并且由于外闲子长度不同,isoform的FPKM也不同,也就是说一个基因可以对应多个FPKM,那这个基因的FPKM应该如何计算?

  33. Reply vidya 4月 11,2017 6:35 上午

    博主您好!打扰了。我用cuffcompare得到combined.gtf文件,从中提出了所需的转录本,请问接下来怎么用cuffdiff做差异表达分析呢?

  34. Reply yy 4月 12,2017 8:41 上午

    博主您好~看了您的博客学习了很多~现在我用cufflinks的流程想要用cuffcompare,将mapping得到的注释信息与参考注释进行比较,找出参考注释信息中没有提及且至少有一个样本的表达值FPKM>=10的转录区域作为novel genes。除了得到的cuffcmp前缀的文件外,还得到了.tmap和.refmap文件,请问我要去哪里寻找novel genes呢?如果对novel genes继续后续做差异性表达,是直接将提出来的novel genes 跑cuffdiff吗?
    实验室里没有做这个部分的师兄师姐,还请您不吝赐教。
    谢谢!

    • Reply admin 4月 13,2017 1:56 下午

      我觉着是这样,首先,我不太同意使用cuffcompare而不用cuffmerge。要知道,cuffmerge是会run Reference Annotation Based Transcript assembler (RABT), 而coffcompare不会。所以在一般情况下,cuffcompare没有cuffmerge准确。二,你可以在输出文件中找到Class code,在那里,你可以找到与reference不同的基因。对于class code, 你可以参考http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/.

      对于差异表达分析,应该把所有的cuffmerge的结果都传进去,否则会影响FDR.

Leave a Reply

  

  

  

%d 博主赞过: