原本,我并無寫這一稿件的想法。主要原因有二:
從事四川樂山服務器托管,服務器租用,云主機,網(wǎng)站空間,域名注冊,CDN,網(wǎng)絡代維等服務。
如果要找合理解釋,那么針對第一點,就是每天仍然有大量新接觸生信數(shù)據(jù)分析的朋友;針對第二點,......在前兩天我推的文稿《零基礎(chǔ)快速完成基因功能注釋 / GO / KEGG / PFAM...》中,評論區(qū)答應了下,閱讀過5000,那就寫一寫富集分析。于是,如果不寫,總是不對。如果要寫,只能現(xiàn)在寫。畢竟有些事情,現(xiàn)在不做,以后真的不會做。
對于這一塊,完全陌生的朋友,尤其是不少生物學背景朋友,有必要溫習一下數(shù)理統(tǒng)計基礎(chǔ)。這一稿件只做原理最簡單的但使用最廣泛其速度最快的Over-Represence Analysis模式的富集分析講演。其他模式,不涉及。
回到主題,先舉個經(jīng)典的抽球例子:
小紅小綠小藍三個人自稱有超能力,可以用手摸摸球就分辨出黑球白球,于是我們找來黑袋子,放100個球,其中20個白球80個黑球,讓三人分別無放回地抽取。
小紅隨機抽出來10個球,其中2個白球8個黑球,情況即,
抽球中白球比例與背景白球比例完全一致,說明小紅抽球結(jié)果隨機。
球放回去,小綠來抽球,抽出來的10個球,其中3個白球7個黑球,情況即,
這是經(jīng)典的抽球案例,抽取到的白球個數(shù)的概率分布為超幾何分布?;诖耍覀兛梢院唵斡嬎愠槿〉奖刃【G抽取到球個數(shù)(或更多即更極端)的概率如何,在 R語言中計算,即
而對于小藍的情況,那么概率如何?
在 TBtools 中也可以計算,只是寫法有點區(qū)別
可以看到,盡管這只是一次抽球,小綠抽球中白球比例(或更極端情況)出現(xiàn)的概率是31.88%+,還是挺高的,于是我們有較高的把握說,小綠嘛,只是走了狗屎運。相反,小藍抽球中白球比例或更極端情況出現(xiàn)的概率幾乎為 0 ,我們幾乎沒啥把握說,小藍走狗屎運....換句話說,我們有理由相信,或許小藍真有抽白球的超能力.....
說了這么多,那么跟基因集合富集分析有啥關(guān)系?....基因集合功能富集分析。那么我們就需要有一個基因集合(如差異表達基因集合或ChIP-seq的Peaks或GWAS定位的系列區(qū)間),還有一個功能標簽(如 生長素信號轉(zhuǎn)導相關(guān) )。于是黑白球案例可以簡單調(diào)整一下。假定現(xiàn)在這個物種一共有100個基因,其中20個基因與生長素信號轉(zhuǎn)導相關(guān),80個沒有注釋到與生長素信號轉(zhuǎn)導相關(guān)(換句話說,約等于無關(guān)),我們做了對植株做了處理,和CK分別測定轉(zhuǎn)錄表達譜,通過差異表達分析,鑒定到10個差異表達基因,其中2個與生長素信號轉(zhuǎn)導相關(guān),而另外8個則沒注釋到生長素信號轉(zhuǎn)導相關(guān),簡單畫一下,即
好,剩下的兩個就不替換了。整體上,ORA模式的富集分析,本身就是經(jīng)典的抽球案例,感興趣的自行替換就可以了。
基本原理,相信都搞清楚了。不過還是有兩三點需要注意:
具體如何做物種所有基因的背景注釋,請參考前述推文《零基礎(chǔ)快速完成基因功能注釋 / GO / KEGG / PFAM...》。
首先,打開 TBtools GO 富集分析界面
整體如上,一共三個文件:
具體示例如下
點擊 Start ,隨后等待即可。完成時會有彈窗提示。查看輸出文件
(寫到這里,突然覺得這些都沒啥意思,不知為何....就不詳細寫了,大伙自己看看列名,猜猜吧)
很多時候,我們會選擇,篩選第一列,只看 Biological Process。一般這些與我們的生物學認知會貼近一些。
基因集合功能富集分析,是一個常常被談起的話題,甚至近期都有不少新方法或算法被提出。感興趣的朋友可以去了解。這份教程,只與大伙說最簡單,但也是使用最為廣泛的一種富集分析模式。無論是不是 TBtools 用戶,理論上來說,都可以輕松理解并掌握,從原理到實踐。
寫到一半,其實我已經(jīng)不想寫了。原因非常簡單,這也是為什么在我之前,并沒有一個人寫出來 TBtools 類似的工具。不是寫不了,而是不想寫。有時候,隨著能力增長和知識積累,往往不再愿意做一些簡單的事情?;蛟S這還涉及到年齡的增長,角色的轉(zhuǎn)變,責任的變化....云云。
小時候,我以為寫 TBtools 玩玩;
后來,我以為我會一直寫下去;
現(xiàn)在,,,,,,
前面我給大家詳細介紹過
?GO簡介及GO富集結(jié)果解讀
?四種GO富集柱形圖、氣泡圖解讀
?GO富集分析四種風格展示結(jié)果—柱形圖,氣泡圖
?KEGG富集分析—柱形圖,氣泡圖,通路圖
? DAVID GO和KEGG富集分析及結(jié)果可視化
也用視頻給大家介紹過
? GO和KEGG富集分析視頻講解
最近有粉絲反映說,利用clusterProfiler這個包繪制GO富集分析氣泡圖和柱形圖的時候,發(fā)現(xiàn)GO條目的名字都重疊在一起了。
氣泡圖
柱形圖
這個圖別說美觀了,簡直不忍直視。經(jīng)過我的認真研究,發(fā)現(xiàn)跟R版本有關(guān)。前面我給大家展示的基本都是R 3.6.3做出來的圖。很多粉絲可能用的都是最新版本的R 4.1.2。
我們知道R的版本在不停的更新,相應的R包也在不停的更新。我把繪制氣泡圖和柱形圖相關(guān)的函數(shù)拿出來認真的研究了一下,終于發(fā)現(xiàn)的癥結(jié)所在。
dotplot這個函數(shù),多了個 label_format 參數(shù)
我們來看看這個參數(shù)究竟是干什么用的,看看參數(shù)說明
label_format :
a numeric value sets wrap length, alternatively a custom function to format axis labels. by default wraps names longer that 30 characters
原來這個參數(shù)默認值是30,當標簽的長度大于30個字符就會被折疊,用多行來展示。既然問題找到了,我們就來調(diào)節(jié)一下這個參數(shù),把他設(shè)置成100,讓我們的標簽可以一行展示。
是不是還是原來的配方,還是熟悉的味道
同樣的柱形圖,我們也能讓他恢復原來的容貌。
關(guān)于如何使用R做GO和KEGG富集分析,可參考下文
GO和KEGG富集分析視頻講解
查看GOplot內(nèi)示例數(shù)據(jù)的格式,對自己的數(shù)據(jù)做處理
觀察結(jié)論:
觀察自己的兩個數(shù)據(jù)表:
table.legend 設(shè)置為T時會顯示表格
本圖中表格和圖例是出圖后剪切拼合而成,沒有用R中的拼圖包
ID轉(zhuǎn)換用到的是 bitr() 函數(shù),bitr()的使用方法:
org.Hs.eg.db包含有多種gene_name的類型
keytypes() :keytypes(x),查看注釋包中可以使用的類型
columns() :類似于keytypes(),針對org.Hs.eg.db兩個函數(shù)返回值一致
select() :select(x, keys, columns, keytype, ...) eg.
函數(shù)enrichGO()進行GO富集分析,enrichGO()的使用方法:
舉例: