基因家族扩增与收缩分析
基因家族扩张和收缩分析
基因家族的扩张和收缩分析一般会使用orthoMCL进行同源基因识别，然后选择直系同源基因进行物种树构建，最后使用CAFE对聚类结果进行基因家族的扩张和收缩分析
直系同源基因鉴定
直系同源基因鉴定网上一般给出了两个软件，一个是orthoMCL，一个是orthofinder。orthoMCL虽然很长时间没有进行过维护更新了，但大家进行基因家族扩张和收缩分析是依然经常性的使用，而orthofinder是16年出现的新软件，本身使用和安装起来更加方便，我也是比较推荐这个
提取最长转录本
进行mcl聚类之前，首先需要挑选每个基因最长的转录本形成一个fa文件，fa文件中是氨基酸序列，不是碱基序列，这里强调一下。然后cafe官方的例子中提供了一个脚本cafetutorial_longest_iso.py，但是这个脚本比较针对他提供的范例fa文件，对于正常分析时的fa文件是不适用的。我在网上找过很多的提取最长转录本的脚本，但是都不具有广泛的适用性，最主要的原因还是几个数据库间对于基因组的注释文件的规则有些不同，而且对于自己组装形成的基因组也有不同，一般自己组装的基因组gff文件是由EVM生成的。而有的生成出来一个基因具有多个转录本，而有的基因组注释并不完全，本身有scaffold构成，一个mrna就是一个转录本，没有重复。所以在提取最长转录本时需要根据实际情况进行提取。这里我的基因组来源有三个，一个是ncbi，一个是ensembl，一个是evm。所以脚本内写了三个模式，不同的gff使用不同的模式。主要思路是对gff进行解析，判断其基因对应的所有mrna，以及所有mrna对应的cds。然后对所有cds进行长度计算，循环遍历每个基因的所有转录本并进行长度比较，保留最长的那条写入pep文件。
需要观察注释的就是解析gff文件这一步，比如对于ensembl注释而言
image.png
判断到mrna之后可以获取到其parent所属的gene是什么，然就从下面的cds可以获得这个cds对应的parent的transcript是什么。而且cds的id就是蛋白id，所以可以使用蛋白序列进行提取，更加方便。
而对于ncbi注释而言，道理是相同的
image.png
不过由于本身id前有部分前缀需要去掉，所以处理时记得注意一下
脚本参考本人github的extract_L_protein.py文件。
同时有时候部分注释会出现点问题，到时候还是要具体情况具体分析，一切根据注释来。然后输入的fa文件记得对蛋白id进行处理，去掉多余的description部分，一切与注释中的结果对应即可
orthofinder运行
orthofinder的运行比较简单。安装好后运行
orthofinder -M msa -T fasttree -f dataset
orthofinder参数并不复杂
-M 指定基因树推断方法，默认‘dendroblast’，可换成‘msa’
-T 物种树推断方法，默认‘fasttree’，可选fasttree, raxml, raxml-ng, iqtree。
-f 表示fa文件夹路径
另外还有一些参数大家可以自己去看--help，比较简单。
fa文件夹就是之前提取的最长转录本的fa文件的文件夹，所有物种放在一起。同时记得将基因id前加上物种名，如human|ENSP00000473549。这样，帮助聚类之后分别基因来源。如何添加更改物种名这里不做赘述，较为简单，用python进行修改即可。
另外我测试时因为未知原因无法使用raxml参数进行建树，使用默认fasttree才可行。
最后结果会在dataset文件夹中生成一个命名类似Results_Jan15的文件夹。里面结果比较多。具体内容参看orthofinder官方文档，或者参考另一篇orthofinder文章。
构建时间尺度物种树
在orthofinder构建成功物种树之后，还不能直接进入cafe中使用，需要使用r8s估计时间尺度。
r8s使用前需要编辑一个input文件，内容和方法参考陈连福老师博客
这里使用cafe在github的脚本cafetutorial_prep_r8s.py。可以直接生成一个input文件
写法如下
python2 ~/tools/CAFE/python_scripts/cafetutorial_prep_r8s.py -i ~/dataset/OrthoFinder/Results_Jan15/Species_Tree/SpeciesTree_rooted.txt -o r8s_ctl_file.txt -s 2335369 -p 'goat,cow' -c 25
-i 表示输入的物种树，我们直接使用生成的物种树，路径如上所示。然后orthofinder的树节点间有1代表分叉节点如((cow:xx.000000,goat:xx.000000)1:x.xxxx，需要删掉1，最后应该是((cow:xx.000000,goat:xx.000000):x.xxxx
-o 表示生成的input文件
-s 表示构建物种树时所用的碱基数量，这个由于是orthofinder构建好了，所以我们需要自己去看其进行比对的序列文件路径为OrthoFinder/Results_Jan15/MultipleSequenceAlignments/SpeciesTreeAlignment.fa，里面包括了所用的物种比对序列，计算其中一个物种的序列长度即可。其实这里部分有点猜想的意思在里面，因为文档就说是进行序列比对的碱基数，我不清楚是所有物种全加上还是其中一个物种的长度。但是观察过范例的碱基数后觉得应该是一个物种，如果是所有的话数字会比较大。
-p 制定某一个分叉节点的时间，如这里制定了羊和牛的分叉时间为25，该数字可以到Timetree进行查询。
然后生成r8s_ctl_file.txt之后如果不满意自动生成的input，还能进行手动修改，比如添加更多分叉节点使结果更加精确，以及分叉节点进行重新命名等，具体参照上面推荐的博客进行修改
准备好input文件后运行r8s，如何安装r8s这里不进行赘述，不难。
~/tools/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt
生成文件r8s_tmp.txt，最后为生成的物种树，设定分离节点会变成字符，这里需要手动将字符删掉。
比如之前设置的分离节点((cow:xx.000000,goat:xx.000000)oatcow:x.xxxx，会多出oatcow，需要将oatcow删掉变成((cow:xx.000000,goat:xx.000000):x.xxxx。
CAFE分析基因家族扩张和收缩
接下来使用CAFE分析基因家族的扩张和收缩
准备cafe.data文件
orthofinder的结果文件OrthoFinder/Results_Jan15/Orthogroups/Orthogroups.GeneCount.tsv就是cafe.data，但是需要稍微修改一下，讲最后一列的total去掉，同时第一列添加desc，内容可以全为null
，cafe的格式要求就是一列desc，一列id，然后物种在该基因家族中的拷贝数。单拷贝基因家族就是该基因家族的所有物种拷贝数为1，之前挑选最长转录本就是因为如果不适用最长转录本会导致大部分基因的拷贝数不为1，使单拷贝基因家族数目不准。修改方法如下
##去掉total
awk 'OFS="\t" {$NF="" ;print $0}' Orthogroups.GeneCount.tsv > cafe.data.1
##增加desc列
sed 's/^/(null)\t&/g' "/home/mwshi/project/xunlu/mcl/compliantFasta/OrthoFinder/Results_Jan15/Orthogroups/cafe.data.1">/home/mwshi/project/xunlu/mcl/compliantFasta/OrthoFinder/Results_Jan15/Orthogroups/cafe.data
##手动把第一列的null改成desc，虽然不改也没关系。不管用什么方法，只要使格式符合cafe输入要求即可
准备cafetutorial_run.sh文件
接下来编辑cafe的运行文件cafetutorial_run.sh，内容如下
#!cafe
load -i cafe.data -t 4 -l log_run1.txt -p 0.05
tree (((chimp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93)
lambda -s -t (((1,1)1,(2,2)2)2,2)
report resultfile
load 以指定要分析的基因家族文件
tree 以指定系统发育树的结构，使用之前生成的时间树即可，记得去掉中间分支节点的字符，否则报错。
lambda 设定不同分支的基因家族进化速率，相同数字表示演化速率相同，不同数字表示演化速率不同，这里感觉数字并不代表速率本身，只是一个编号，我自己运行时由于不知道互相演化速率差异所以全部填1进行运行，软件自身也会对演化速率进行估计。
report 表示返回的结果文件
编辑好后运行cafe
cafe cafetutorial_run1.sh
结果储存在resultfile.cafe中。
cafe结果可视化
这一步对我非常困难，因为在网上搜cafe可视化发现并没有什么资料，都是说自己进行结果提取和画图。
这里有两种方法，一是使用cafe自带脚本计算出扩张和收缩数目后自行绘图，比如ggtree或者itol手动绘制等。
方法如下
python2 ~/tools/CAFE/python_scripts/cafetutorial_report_analysis.py -i resultfile.cafe -o summary_run1
python ~/tools/CAFE/python_scripts/cafetutorial_draw_tree.py -i summary_run1_node.txt -t '(((chimp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93) ' -d '(((chimp<3>,human<3>)<3>,(mouse<3>,rat<3>)<3>)<3>,dog<3>) ' -o summary_run1_tree_rapid.png
其中cafetutorial_report_analysis.py生成结果文件中summary_run1_node.txt中包含每个节点扩张和收缩的数目。cafetutorial_draw_tree.py简单的对结果进行绘图，以上参数在resultfile.cafe中找。默认绘制扩张基因数目，可以添加-y Contractions来绘制收缩基因数目，然后根据这些使用其他工具绘制更美观的图片
第二种方法是我找到的唯一一个工具直接对cafe结果进行可视化的工具，CAFE_fig。
安装时其主页现实安装ete3 3.0.0b35版本，不用管它，安装最新版本。他提示安装低版本是英文有部分人会因为其ete3无法使用PyQt5的内容所以必须要降级。但是如果选择使用3.0.0b35版本代表必须将PyQt5降级到PyQt4，也比较麻烦。而且我自己测试时全部使用最新版并没有出现问题，所以安装时不需要降级。
python3 ~/tools/CAFE_fig/CAFE_fig.py resultfile.cafe -pb 0.05 -pf 0.05 --dump test/ -g pdf --count_all_expansions
直接运行即可，另外，由于本人物种树时间长度不足，导致画出的图非常小。这里提供两种方法。第一种是对图片的比例尺进行估计，比如想让树边长5倍，则在cafe时对树节点值*5，扩大长度，然后画图完成后使用ai对比例尺的值/5修正即可。第二种方法是我阅读其代码，发现他使用的ETE3进行绘图，参照ETE3说明文档之后只需要在CAFE_fig.py中对pixels_per_mya进行修正。源码中是1.0，如果想让长度变成5倍则更正为5.0即可
image.png
然后由于未知原因，我画出的图片没有颜色，暂时还不清楚为什么。
以上就是我弄出的基因家族扩张和收缩分析流程，因为本身在这一块还属于新手上路，可能整个流程还有很多缺陷需要弥补，希望各位大佬不吝赐教。然后我也不清楚是否还有更简单，更高效的方法流程，如果有大佬愿意指点指点，感激不尽。
