最近剛返校溯饵，事情比較多七冲，每天也很忙，之前寫(xiě)的《基因家族擴(kuò)張與收縮分析及物種進(jìn)化樹(shù)構(gòu)建（上）》也一直沒(méi)來(lái)得及更新殿较，缺少cafe輸出結(jié)果的解讀及后面的可視化耸峭。最近在簡(jiǎn)書(shū)上也收到了不少私信，也很高興自己利用閑暇時(shí)間寫(xiě)的教程能對(duì)大家的生信分析起到略微的幫助淋纲，這也是我開(kāi)通簡(jiǎn)書(shū)的初衷劳闹。我自己其實(shí)也不是生信專(zhuān)業(yè)出身，也是由于課題需要才接觸的帚戳，從實(shí)驗(yàn)室的服務(wù)器環(huán)境配置玷或、軟件安裝與調(diào)試，到后面接觸項(xiàng)目片任，自己做分析偏友，到現(xiàn)在也算是入門(mén)了吧。這一路走來(lái)对供，遇到的坑不少位他，踩的坑也挺多，也是在網(wǎng)上看其他人寫(xiě)的教程产场，一步步嘗試鹅髓，一次次解決報(bào)錯(cuò)，這樣摸爬滾打走過(guò)來(lái)的京景，單純的野路子窿冯。一方面也是受此影響，另一方面也是為了總結(jié)和整理自己的學(xué)習(xí)筆記确徙，所以寫(xiě)了這些文章醒串。由于寫(xiě)每一篇文章执桌，都是用的真實(shí)數(shù)據(jù)，其中這篇《基因家族擴(kuò)張與收縮分析及物種進(jìn)化樹(shù)構(gòu)建》也是選的10個(gè)園藝果樹(shù)物種的基因組外加擬南芥和水稻芜赌，從頭邊做分析仰挣，邊寫(xiě)的，也是確保每一步的代碼運(yùn)行的穩(wěn)定且可靠缠沈，而且每一步的報(bào)錯(cuò)也能詳細(xì)的指出膘壶。其中有一步是物種進(jìn)化樹(shù)的構(gòu)建，我采用的是iQtree做的洲愤，利用的是最大似然法構(gòu)建的STAG有根樹(shù)颓芭，所以運(yùn)行較慢，這一步足足跑了兩個(gè)星期禽篱！

最終結(jié)果如下：(Oryza_sativa:0.381481,(Diospyros_oleifera:0.269879,(Vitis_vinifera:0.190065,((Prunus_persica:0.0744669,(Malus_domestica:0.0156371,Pyrus_communis:0.0223021):1):1,((Dimocarpus_longan:0.177168,((Citrus_clementina:0.00928013,Citrus_sinensis:0.00417004):1,Citrus_grandis:0.00781403):1):1,(Arabidopsis_thaliana:0.44985,Durio_zibethinus:0.205998):1):1):1):1):1)

12-spp-STAG-tree

得到了物種進(jìn)化樹(shù)之后畜伐，這里用R8S提取時(shí)間樹(shù)，此步參照上一篇躺率；直接放結(jié)果：
((((((Malus_domestica:0.722166,Pyrus_communis:0.722166):28.002356,Prunus_persica:28.724522):46.395611,((((Citrus_clementina:0.154051,Citrus_sinensis:0.154051):14.174580,Citrus_grandis:14.328630):18.378373,Dimocarpus_longan:32.707003):23.254341,(Arabidopsis_thaliana:17.372300,Durio_zibethinus:17.372300):38.589044):19.158789):17.535799,Vitis_vinifera:92.655932):23.603028,Diospyros_oleifera:116.258960):35.741040,Oryza_sativa:152.000000)

有了這兩個(gè)結(jié)果后就可以配置CAFE的輸入文件了

cafetutorial_run.sh

接下里運(yùn)行cafe，進(jìn)行基因家族擴(kuò)張與收縮分析：
nohup cafe cafetutorial_run.sh 2> cafe.log &

得到的cafe輸出結(jié)果如下：

out.cafe

官方文檔-cafe輸出文件格式說(shuō)明

cafe輸出結(jié)果

根據(jù)官方文檔介紹万矾，物種名后下劃線“_”后面的數(shù)字即為該物種在該基因家族的基因數(shù)目悼吱，那么收縮和擴(kuò)張時(shí)如何來(lái)計(jì)算和判定的呢？先看第四行的結(jié)果：

cafe輸出結(jié)果第四行

例如：我們想知道甜橙中基因家族擴(kuò)張和收縮的情況良狈，我們注意到甜橙所在的節(jié)點(diǎn)是<8>后添，在OG0000001中，它的基因數(shù)目是39薪丁，而與它相近的支節(jié)點(diǎn)為<7>遇西，其下劃線“_”對(duì)應(yīng)的數(shù)字為0，如果<8> - <7> 結(jié)果是大于零的严嗜，那么針對(duì)甜橙而言粱檀，該家族就是擴(kuò)張的，反之如果小于零漫玄，那么就是收縮的茄蚯，如果等于零，那就是“no change”睦优。針對(duì)此渗常，編寫(xiě)腳本很容易提取相關(guān)的信息。

最后就是可視化繪圖了汗盘，這里我們用“CAFE_fig”來(lái)化皱碘，cafe的輸出結(jié)果可直接做輸入文件：

python3 CAFE_fig.py out.cafe -pb 0.05 -pf 0.05 --dump test/ -g svg --count_all_expansions

這里的輸出文件為svg格式，可直接導(dǎo)入到AI里進(jìn)行編輯和美化隐孽。最終呈現(xiàn)結(jié)果如下：

12-spp-STAG-tree

以上