最近剛返校，事情比較多，每天也很忙，之前寫的《基因家族擴(kuò)張與收縮分析及物種進(jìn)化樹(shù)構(gòu)建（上）》也一直沒(méi)來(lái)得及更新，缺少cafe輸出結(jié)果的解讀及后面的可視化。最近在簡(jiǎn)書上也收到了不少私信，也很高興自己利用閑暇時(shí)間寫的教程能對(duì)大家的生信分析起到略微的幫助，這也是我開(kāi)通簡(jiǎn)書的初衷。我自己其實(shí)也不是生信專業(yè)出身，也是由于課題需要才接觸的，從實(shí)驗(yàn)室的服務(wù)器環(huán)境配置、軟件安裝與調(diào)試，到后面接觸項(xiàng)目，自己做分析，到現(xiàn)在也算是入門了吧。這一路走來(lái)，遇到的坑不少，踩的坑也挺多，也是在網(wǎng)上看其他人寫的教程，一步步嘗試，一次次解決報(bào)錯(cuò)，這樣摸爬滾打走過(guò)來(lái)的，單純的野路子。一方面也是受此影響，另一方面也是為了總結(jié)和整理自己的學(xué)習(xí)筆記，所以寫了這些文章。由于寫每一篇文章，都是用的真實(shí)數(shù)據(jù)，其中這篇《基因家族擴(kuò)張與收縮分析及物種進(jìn)化樹(shù)構(gòu)建》也是選的10個(gè)園藝果樹(shù)物種的基因組外加擬南芥和水稻，從頭邊做分析，邊寫的，也是確保每一步的代碼運(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)
得到了物種進(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的輸入文件了
接下里運(yùn)行cafe，進(jìn)行基因家族擴(kuò)張與收縮分析：
nohup cafe cafetutorial_run.sh 2> cafe.log &
得到的cafe輸出結(jié)果如下：
根據(jù)官方文檔介紹，物種名后下劃線“_”后面的數(shù)字即為該物種在該基因家族的基因數(shù)目，那么收縮和擴(kuò)張時(shí)如何來(lái)計(jì)算和判定的呢？先看第四行的結(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ì)此，編寫腳本很容易提取相關(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é)果如下：
以上
