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研究會,獲取知識,開拓事業。