最近刚返校，事情比较多，每天也很忙，之前写的《基因家族扩张与收缩分析及物种进化树构建（上）》也一直没来得及更新，缺少cafe输出结果的解读及后面的可视化。最近在简书上也收到了不少私信，也很高兴自己利用闲暇时间写的教程能对大家的生信分析起到略微的帮助，这也是我开通简书的初衷。我自己其实也不是生信专业出身，也是由于课题需要才接触的，从实验室的服务器环境配置、软件安装与调试，到后面接触项目，自己做分析，到现在也算是入门了吧。这一路走来，遇到的坑不少，踩的坑也挺多，也是在网上看其他人写的教程，一步步尝试，一次次解决报错，这样摸爬滚打走过来的，单纯的野路子。一方面也是受此影响，另一方面也是为了总结和整理自己的学习笔记，所以写了这些文章。由于写每一篇文章，都是用的真实数据，其中这篇《基因家族扩张与收缩分析及物种进化树构建》也是选的10个园艺果树物种的基因组外加拟南芥和水稻，从头边做分析，边写的，也是确保每一步的代码运行的稳定且可靠，而且每一步的报错也能详细的指出。其中有一步是物种进化树的构建，我采用的是iQtree做的，利用的是最大似然法构建的STAG有根树，所以运行较慢，这一步足足跑了两个星期！
最终结果如下：(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
得到了物种进化树之后，这里用R8S提取时间树，此步参照上一篇；直接放结果：
((((((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)
有了这两个结果后就可以配置CAFE的输入文件了
cafetutorial_run.sh
接下里运行cafe，进行基因家族扩张与收缩分析：
nohup cafe cafetutorial_run.sh 2> cafe.log &
得到的cafe输出结果如下：
out.cafe
官方文档-cafe输出文件格式说明
cafe输出结果
根据官方文档介绍，物种名后下划线“_”后面的数字即为该物种在该基因家族的基因数目，那么收缩和扩张时如何来计算和判定的呢？先看第四行的结果：
cafe输出结果第四行
例如：我们想知道甜橙中基因家族扩张和收缩的情况，我们注意到甜橙所在的节点是<8>，在OG0000001中，它的基因数目是39，而与它相近的支节点为<7>，其下划线“_”对应的数字为0，如果<8> - <7> 结果是大于零的，那么针对甜橙而言，该家族就是扩张的，反之如果小于零，那么就是收缩的，如果等于零，那就是“no change”。针对此，编写脚本很容易提取相关的信息。
最后就是可视化绘图了，这里我们用“CAFE_fig”来化，cafe的输出结果可直接做输入文件：
python3 CAFE_fig.py out.cafe -pb 0.05 -pf 0.05 --dump test/ -g svg --count_all_expansions
这里的输出文件为svg格式，可直接导入到AI里进行编辑和美化。最终呈现结果如下：
12-spp-STAG-tree
以上
网友评论