Orthofinder+MCMCTree+Cafe完成系统发育及基因家族扩张收缩分析
一、数据获取
物种间的基因家族扩张收缩分析一般以蛋白序列来做。目标物种的蛋白序列可以从ncbi、esemble、JGI等数据库获取。
二、数据处理
各个数据库间的注释信息等可能存在差异，获取后需要提取最长蛋白序列，简化id，否则会影响提取同源单拷贝基因等后续分析。关于提取最长蛋白序列，网上脚本很多。这个步骤非常重要，为了方便后面分析需要处理好。
三、Orthofinder同源物分析
运行完后，提取对应的单拷贝序列，使用mafft完成多序列比对，然后将所有的序列串联在一起，使用Raxml-HPCPTHREADS建树，获取同源单拷贝构建的最大似然树。当然以上操作也可以通过设置Orthofinder一步到位。例如添加使用-M msa - T raxml-NG等参数。
四、MCMCTree构建时间分化树
该步骤利用构建的物种树文件，通过添加时间分化定点（一般最好大的分支上都要有定点），对整个物种树的时间分化进行推断。http://www.timetree.org/可以查询两两物种的分化时间。mcmctree流程参照官方文档http://abacus.gene.ucl.ac.uk/software/MCMCtree.Tutorials.pdf即可。
五、CAFE基因家族扩张收缩分析
经过Orthofinder过后，我们得到了所有物种同源的基因家族，应用CAFE可以对基因家族的历史进行推断。CAFE应用人口出生死亡模型对每一个基因家族的历史进行推断，判断其经历扩张收缩还是保持不变。经过这一分析我们可以得知我们关注的物种中哪些基因家族发生了扩张和收缩，这些基因家族可以进行进一步的分析判断与哪些生物学过程相关。CAFE的使用参照官方文档https://hahnlab.github.io/CAFE/src_docs/html/index.html 以及https://github.com/hahnlab/cafe_tutorial/blob/main/cafe_tutorial.pdf需要注意的是MCMCTree进行时间推断时设置的分化时间是以100个MYA为单位的，而在输入到CAFE中的时间分化树需要改为以MYA为单位。
另外尤其需要注意一点，CAFE的关键参数有一lambda，输入文件有二：其一是时间分化树，上面已经说明了注意事项；其二是GeneCounts.tsv文件，做基因家族扩张收缩分析应该将该文件N0节点上不存在的OG过滤掉（也就是根节点排序过滤掉为0的那些OG）。这一点CAFE4似乎不会替我们完成，输出结果和过滤后并不一致，CAFE5应该是会自动过滤的。
六、结果可视化
CAFE结果可视化使用https://github.com/LKremer/CAFE_fig 缺什么包安装什么包即可，需要注意的是如果确认了包都安装成功后，仍出现找不到包的报错时，需要确认一下包安装在了什么位置，使用对应版本的python的绝对路径运行脚本即可。如~/biosoft/envs/python3.8/bin/python3 CAFE_fig.py。或者可以在安装时使用绝对路径的pip3去安装对应的python包，最后调用对应位置的python程序就不会出现找不到包的报错了。其他细节可翻阅其他相关文档。
七、总结
CAFE5更好用，命令可以更简单，不需要再准备lambda树等麻烦操作。一行就可以搞定。
#Orthogroups.GeneCounts.tsv是Orthofinder的输出结果，需提取处理格式即可。
#Mcmctree.out.dated.tre是mcmctree的结果，提取出树内容保存即可。
~/biosoft/CAFE5/bin/cafe5 -i Orthogroups.GeneCounts.tsv -t Mcmctree.out.dated.tre