如何計算cuffdiff中的FPKM值(含Q&A)

2021-01-15 生物狗課堂

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的公式就可以從字面上寫成:


,其中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中是這樣描述的:


其中beta和gamma都是likelihood function中的參數。lt被定義為


因為,假設transcript t的長度為l(t),那麼隨機分布時在某一位置t出現長度為k的片段的概率為


在cuffdiff中,它會將同一組中所有的樣品試為同一來源樣品處理,這就是為什麼同一組內如果有三個樣品的話,最終得到FPKM值並不是三個樣品單獨運算得到的FPKM值的平均值。


Q&A


Q:「在cuffdiff中,它會將同一組中所有的樣品試為同一來源樣品處理,這就是為什麼同一組內如果有三個樣品的話,最終得到FPKM值並不是三個樣品單獨運算得到的FPKM值的平均值。」
您好,我曾經試著分開運算每個樣品,確實發現總的結果不是平均值。可是如果將所有樣品作為一個sample來處理,那後面的p值是怎麼得到的呢?我的數據裡有些基因在一個sample裡表達量是0,但在同一組另一個sample裡FPKM卻有幾百。不知道這個是否正常?


A:你問到有些基因在一個sample中的表達量與其它同組的sample的表達量差別很大,是否正常。這完全有可能是正常的。我時常遇到這樣的數據,但最後得到的結果驗證都還不錯。這就是生物樣品的複雜性,即使平行實驗,也可能差異很大。


Q:就是在運行玩cuffdiff以後發現對應FPKM最高值的大部分基因居然都沒有gene name。我的預期是高表達量的區域應該對應著house keeping gene,那麼怎麼會還沒有name呢?另外就是沒有gene name的轉錄區域是不是就是還未發現的新的轉錄區呢?謝謝!

A:這個問題比較複雜。一個是檢查你的注釋文件是否正確,是否在cuffdiff時正確使用。第二就是查看那些沒有gene name的轉錄區域是否真的沒有已知注釋。但是沒有gene name的轉錄區並不意味著一定就是新的未知轉錄區。


Q:關於注釋文件是否正確,確實,在我run cuffmerge和cuffdiff的時候都出現了「warning: Couldn’t find fasta record for 『chr13_random』」…還有好幾個chromosome(但不是所有)也有一樣的問題。我試著從不同的地方下reference,但結果都是這樣。我在網上看到有人和我遇到了相似的問題,但都沒有好辦法解決。不知道這個問題會不會有影響。

A:這些警告信息都是可以接受的。因為在那麼random chromosome上本來就沒有什麼好的注釋。


Q:那麼關於我的另一個問題,就是既然cuffdiff裡的value具體是怎麼算的呢?當我對每個sample單獨運行cuffdiff時,發現個別基因在整個運行時得出的value比每一個sample單獨運行cuffdiff時的value都要小,說明總的value確實不是平均值。 還有那個p值是基於每個sample的FPKM算出來的嗎?

A:在其說明中這樣寫的: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.


Q:如何能知道cuffdiff使用了多少tophat結果中的reads呢?也就是說,比對的結果中有多少reads對於計算差異表達是useful的呢?也就是說,如何count有多少reads能比對到cuffdiff所使用的GTF上?


A:我想你可以在read_groups.info以及genes.read_group_tracking文件中找到自己需要的答案。


Q:我發現了一個現象,對於A、B兩個樣品,針對特定一個GTF,有時候,A比對上的fragment比B的多(多1%),但分別對A和B的FPKM求和後,居然發現A的FPKM之和小於B(少了近20%),不知您怎麼看這種現象?


A:有很多種原因可能導致這一結果。1.非unique mapped reads。2.single ends數多。3.cufflinks發現了新的transcripts。我想到的暫時就這麼多了。


Q:就是cuffdiff的作用到底是幹嘛的,它不是進行sam和sam之間的比對嗎,為什麼在輸入文件中還要提供一個.gtf文件?這個.gtf文件起什麼作用呢?
我是個小白,還希望您不吝賜教。謝謝。


A:對於map好的測序結果,一般只包括reads在基因組中的位置,並沒有注釋到每個基因有多少個reads。在比較時,輸入gtf文件讓cuffdiff可以依據一定的模型來計算對應到每個具體的基因上會有多少個reads。



Q:您好,我想問下,可以直接從merged.gtf文件中利用什麼工具計算出FPKM的值嗎?或者是,前面使用cufflink後的結果文件transcript.gtf中是有FPKM的,但這要怎麼一一對應到merged.gtf文件中去呢?其實我現在思路也有點亂,不知道該怎麼問比較合適,或許您能理解我的意思,麻煩了!謝謝~

A:你可以使用htseq-count來獲得一個計數表。如果需要FPKM值的話,你可以在cuffdiff的結果中找到,有一些track文件,其中有每個樣品的FPKM值。


Q:我想問下,用cuffdiff分析差異表達基因時,如果有4個樣本,兩個病例,兩個對照,分析時的merged.dtf是4個樣本合在一起的轉錄本麼?要比較病例和對照之間的差異表達,cuffdiff語句後面跟的是4個樣本的.bam文件麼?

A:看你的需要。如果你需要發現新的轉錄本的話,是需要4個樣本合在一起的轉錄本。如果是其它情況,有可能只需要已知的轉錄本。在比較時,應該把所有涉及到的bam文件都傳進去。


Q:4個樣本中2個對照,兩個病例,要找出他們的差異表達基因,但cuffdiff運行得出來的結果gene_exp.diff文件中是分別作了4個樣本間兩兩的差異表達,怎樣可以得到病例組和對照組的差異表達基因啊

A:你注意,將bam文件傳遞給cuffdiff時,相同組使用逗號(,)連接即可。你可以仔細閱讀一下說明文檔。


Q:這樣做出來的差異表達基因以Q值<0.05和log2(fold change)絕對值大於2刪選,結果用做通路富集,可是沒有富集的通路是怎麼回事啊

A:這很可能是正常的。你現在可以試試GSEA分析。GSEA就是為這種情況而生的。


Q:我想問個問題。cufflinks在計算FPKM的時候,如果read1比對上chr1,read2比對上chr9,這樣這麼算呢?
看了比對結果統計,感覺Reads mapped in proper pairs不是很高,這麼多這樣的情況cufflinks是怎麼算的啊?

A:paired end在map的時候就會要求輸入兩個reads之間可能的距離,所以並不會出現你所謂的read1在一個chromosome,而reads2在另一個chromosome。


Q:請問一下鏈特異性建庫怎麼找反義轉錄本。cufflinks可以做嗎?

A:如果是特異性建庫的話,是可以的。但到了後期注釋,還是需要自己手工檢查。


Q:您好!我是用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的效果?

A:對於樣品特別多,請使用cluster或者雲計算。
1. cufflinks在設計之初需要考慮兩種情況,一種是發現新的轉錄本,一種是不需要發現新的轉錄本。如果直接是第二種,可能就不會經過cufflinks以及cuffmerge,所以在參數設置時,不可能假設你已經有了這些輸入。第二,cufflinks給出的結果沒有經過cuffmerge這一步,不具有比較的意義。
2. 正確。但是,FPKM本身是會依據cuffmerge的結果而變化的。有時候因為比對的樣品不同,對於相同的樣品FPKM值也會不同(當然差別不是很大),你做得多了會慢慢體會。所以,你可以這麼理解,但又不全然如此。
3. 你可以查看htseq-count + R/bioconductor:::edgeR(or DESeq)。但是沒有cluster的幫助,你計算量很大的情況下,一定非常慢。即使有cluster,你也要把心態放寬,1,2個星期是正常的,我時常要等1個月才能看到結果。


Q:假如我不考慮新的轉錄本。可以用cufflinks 的FPKM結果,自行做樣本間normalization 和差異表達檢驗。這個與cuffdiff出來的結果差別嗎?
因為HTseq-DESeq這個流程只能給出差異表達基因,而cuffdiff還能給出差異表達的isoform,所以我特別想利用cufflinks isoform.tracks 這個結果,但是沒有計算資源支持我用cuffdiff。

A:你可以考慮使用DEXseq。這個可以得到差異表達的isoform,不過可能會比較複雜。而且在性能上也不一定會比cuffdiff強。還有,如果使用DESeq,還是使用raw reads count會好些。這些信息可以在track file裡得到。


Q:我想問兩個問題:(1)如何用cufflink軟體確定一個新基因。具體的code或是從哪些結果中去分析。
(2)如何單獨計算一個沒有gene_name的基因的FPKM。

A:您如果運行過cufflinks的話,我想您可以很容易從輸出文件夾中找到你所有的答案。如果不能滿足你的需要的話,可以試試cummeRbund。



Q:您好,請問
1、用cuffdiff計算出來的gene_exp.diff中每一個基因的value值是不是就是FPKM?
2、「在cuffdiff中,它會將同一組中所有的樣品試為同一來源樣品處理,這就是為什麼同一組內如果有三個樣品的話,最終得到FPKM值並不是三個樣品單獨運算得到的FPKM值的平均值。」
同一來源樣品處理,是不是就是說把3個重複的合併成一個樣品,合併之後再計算FPKM,這樣的話p-value又該怎麼算?還有一個q-value是什麼呢?
3、雖然我看了您在回答中關於p-value的一段引用,對於它的計算還是覺得挺模糊,請指教。

A: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.


Q:我也是用的tophat-cufflinks-cuffdiff這個流程,請問一下,cuffdiff結果裡一些轉錄本的FPKM值都為0可能是什麼原因?

A:都是0是因為它本身有很少的幾條reads,當計算FPKM值時小數點截取之後相當於0了。


Q:現在有個問題:我用cufflinks 來計算FPKM值,但結果和Raw Count(由htseq-count計算) 出入很大,比如:四個樣本的RawCount分別是0,0,5744,4606 而RPKM卻是0,0,0,402.494。這齣入也太大了啊,請教我該怎麼辦。重複了幾次都是這樣,但也沒見報錯。
還有博主還有沒有其他計算RPKM和FPKM的工具,給推薦一下,非常感謝。

A:其它的工具我也沒有了,主流就是這兩個。RPKM值和rawCount之間的轉換的確挺難的。如果你給的值如果不是0多輸入了一個的話,我覺得是非常正常的結果,沒有出入很大啊。你需要正確理解RPKM,它和raw count不完全一樣,它是標準化之後的值。


Q:我想問一下cuffdiff的結果中fold change為多少的時候可以確定基因在一個樣本中是特異的?還有cuffdiff結果中的P value是怎麼算的呢?泊松分布?另外,如果在多個樣本(比如6個)之間篩選特異的基因該怎麼做呢?是將兩兩比較的結果合併到一起再篩選?這裡面應該有一個normalization的問題吧?十分感謝!

A:p值的計算比你想像的要複雜,不是簡單的泊松分布就可以的。normalization的過程在最初的FPKM值的時候就已經做過了。具體請閱讀原始文獻。


Q:我想請問一下我按照nature protocol上的那篇文章的步驟用tophat-cufflinks-cuffdiff這個流程來找差異表達的基因,但是最後得到的結果gene_exp.diff中significant沒有為yes的,可能有什麼原因呢?
謝謝博主!

A:原因有很多種,最可能出現的問題在於平行樣品之間差異過大,或者比對樣品之間差異過小。其它原因,比如reads數不足,或者說基因組覆蓋度不夠,或者你的腳本中有小錯誤等等。


Q:您說的那些原因應該都不是,因為數據還有實驗步驟我都是按照《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)的呢?謝謝!

A:它先通過 Poisson model來計算出一個pvalue,然後用這個pvalue較正一下得到qvalue,也就是FDR,然後使用fold change 和qvalue來判斷是否significant.


Q:cuffdiff結果如何做樣本相關性重複分析?還有cuffdiff結果不會顯示每一個重複的表達倍數變化,有沒有可以計算每個重複樣本的表達變化?謝謝

A:cuffdiff在這方面並不擅長。但是它還是有輸出每個重複樣品的結果的,你多注意一下那些tracking文件。


Q:關於FPKM,我一直都有個疑問:FPKM計算得到的是每個轉錄本(isoform)的表達量,一個基因可以有多個轉錄本,這個基因的每個轉錄本(isoform)分別有個FPKM,並且由於外閒子長度不同,isoform的FPKM也不同,也就是說一個基因可以對應多個FPKM,那這個基因的FPKM應該如何計算?

A:基因的範圍是指所有isoform的範圍。


本文來自糗世界的博客,謝謝作者辛勤勞動和無私分享。由基因學院(jiyinxueyuan.com)整理。


生物狗課堂,生物人的學習天地。跟大咖一起學習。  加入itbiomed研究會,獲取知識,開拓事業。

相關焦點

  • Git命令解析-patch、apply、diff
    Part 2我們在使用版本控制工具時,總會花費很多時間來處理diff,比如檢查正在進行的未提交的工作,查看單個提交中發生了什麼變更,在執行合併之前比較兩個分支,等等。diff是版本控制的核心概念,但可能大多數使用者沒有考慮過它是如何生成的。請思考一下如何編寫一個函數來計算diff。很容易發現,更準確的需求是只顯示發生了什麼變化,忽略其他保持不變的部分。
  • 電感Q值計算公式及Q值影響因素詳解
    打開APP 電感Q值計算公式及Q值影響因素詳解 佚名 發表於 2016-11-30 10:52:29 對調諧迴路中的電感線圈,Q值要求較高,因為Q值越高,迴路的損耗就越小,迴路的效率就越高;對耦合線圈來說,Q值可以低一些;而對於低頻或高頻扼流圈,則可以不做要求。 實際上,Q值的提高往往受到一些因素的限制,如導線的直流電阻、線圈骨架的介質損耗、鐵心和屏蔽引起的損耗以及高頻工作時的集膚效應等。因此,線圈的Q值不可能做得很高,通常Q值為幾十至一百,最高也只有四五百。
  • 如何利用PyTorch中的Moco-V2減少計算約束
    介紹SimCLR論文(http://cse.iitkgp.ac.in/~arastogi/papers/simclr.pdf)解釋了這個框架如何從更大的模型和更大的批處理中獲益,並且如果有足夠的計算能力,可以產生與監督模型類似的結果。但是這些需求使得框架的計算量相當大。
  • MATLAB符號計算(收藏版)
    MATLAB符號計算1 符號計算基礎1.1符號運算1.2 符號對象1.3 符號表達式中變量的確定2 符號導數及其應用3 符號積分4 級數符號求和>pi2=pi; r1=8; r2=2; r3=3;                         % 定義數值變量sin(pi1/3)                                    % 計算符號表達式值  sin(pi2/3)                                     % 計算數值表達式值
  • 以前老外喜歡把演講詞記在袖口上,那off the cuff是啥意思
    Off the cuff的中文意思是即席的;即興的。據說這個說法起源於20世紀30年代,當時的人有個習慣,他們會把事情記在袖口上,那時的袖口是塑膠做的。例如:Be careful of any off the cuff remarks you make to reporters.
  • RPKM, FPKM, TPM有什麼區別?
    與q-PCR類似,轉錄組數據基因表達量的衡量也是採用相對定量的方法。那麼,到底該如何定量呢?假設通過比對得到了A B C D 四個基因在 3個樣本中的 read counts 如下表,那麼,用比對到參考基因的Reads個數來衡量不同樣本基因的表達量可以嗎?可以認為樣本1中gene B比gene C表達量高嗎?或者gene B在樣本3中比樣本1表達量高嗎?
  • r語言檢驗時間序列平穩性_r中時間序列平穩怎麼檢驗 - CSDN
    我們可以使用R「forecast」包中的「forecast.HoltWinters()」函數對更多時間點進行預測。要使用forecast.HoltWinters()函數,我們首先需要安裝「預測」R包(有關如何安裝R包的說明,請參閱如何安裝R包)。
  • 如何計算玻璃的K值?
    節能玻璃的重要參數是K值,可以通過計算和軟體的方式對玻璃的K值進行計算,沒有必要對玻璃做完後進行檢測,通過軟體和計算的方法可以節省實際檢測的費用並且可以獲得較準確的數據。下面舉例計算玻璃的K值。一、計算使用公式法。
  • 如何計算共射極放大電路的各個參數
    如何計算共射極放大電路的各個參數。很實用。 實際在應用過程中,就是需要確定上述各個電阻的取值。 Icq =( VCC – Vcq)/Rc = 5mA ; Icq≈Ieq = 5mA; Veq = Ieq * Re = 0.5V; 假設三極體的Vbe = 0.7V,所以 Vbq = Vbe + Veq = 1.2V; 注意,在模電書中會有一個條件,若是三極體的輸入阻抗Rin-resistance ≈ hFE *
  • API 測試工具 Hitchhiker v0.7,Schedule 的對比 diff
    支持以diff方式查看Schedule的對比結果Hitchhiker的Schedule是支持不同環境的數據對比的,不過之前只是把兩邊的response和對比結果給出來,想要知道有哪些不同的話還需要藉助其他diff工具來對比,比較麻煩。
  • 角度為自然數度數的sin值是否都能算出精確值?比如sin1°等
    我們很容易計算sin18°的精確值,它可以用含有根號的一個代數式表示。  .  那麼問題來了,還有哪些角度為自然數度數的sin值能算出精確值?  根據三角函數的特性,我們只需要考慮n=1~90自然數時,sin(n°)是否能計算就可以了。  .
  • 如何一鍵計算OR值?
    OR(比值比,odds ratio),其是病例對照研究(回顧性研究)中常用的指標,用於測量暴露因素與疾病因素之間關聯強度;其實際意義為:暴露組疾病危險度 是 對照組疾病危險度的多少倍。危險度=發病數 / 非發病數。當OR值>1時,說明暴露使疾病危險度增高,是疾病的危險因素。
  • 如何屏蔽掉spice model中已經提取的寄生電容?
    在寄生參數提取中,除了需要關注金屬導體的寄生電容外,還需要分析mos器件的節點寄生參數提取是否準確。下面是一個典型的mos管器件:該mos管的4個節點D, G, S, B的寄生電容值是否準確呢?G到B的電容比剛才的值小了,說明在mos管區域的電容被屏蔽了,只剩下Gate露頭的電容了。這個符合預期。但是,S到B的電容和D到B的電容沒有屏蔽掉,該如何設置呢?
  • 生物統計(7)——R中dnorm, pnorm, qnorm與rnorm的區別
    我們知道,在R中,與正態分布(或者說其它分布)有關的函數有四個,分別為dnorm,pnorm,qnorm和rnorm,其中,dnorm表示密度函數,pnorm表示分布函數,qnorm表示分位數函數,rnorm表示生成隨機數的函數。