【比较基因组】基因家族扩张与收缩---CAFE篇

====环境和测试数据准备=====

需要先安装OrthoFinder，这个参考前面一个帖子。然后今天我们主要来测试CAFE。

 安装也比较简单：

下载安装包之后

https://github.com/hahnlab/CAFE/releases/download/v4.2.1/CAFE-4.2.1.tar.gz

cd CAFE

./configure

make

make install prefix=文件夹

注：在release文件夹中成功安装了cafe，可以把它加入自己的环境变量。

测试数据准备：一共下载了mouse, rat, cow, horse, cat, marmoset,macaque, gibbon, baboon, orangutan, chimpanzee, 和 human 12个物种的蛋白数据。（我是另外一个博主那里下载的：从https://share.weiyun.com/5ZIjBg8 (密码：3jzdpm)下载。）

//解压缩一下就可以了

tar xf twelve_spp_proteins.tar.gz

for i in `ls -1twelve_spp_proteins/*.tar.gz`; do tar xf $i -C twelve_spp_proteins; done

rm twelve_spp_proteins/*.tar.gz

===识别基因家族=====

识别物种内和物种间的基因家族分为如下四步：（当然也可以）

 

1、仅保留每个基因中有代表性的转录本，去除可变剪切和冗余基因

2、建立BLAST数据库，使用blastp进行 all-by-all 的比对

3、使用MCL基于blastp结果进行聚类，基因序列相似的通常是一个基因家族

4、解析MCL的输出结果，用作CAFE的输入

第一步：将所有最长的转录本合并成单个文件。提取每个基因中最长的转录本，然后合并成单个文件。

python python_scripts/cafetutorial_longest_iso.py -d twelve_spp_proteins/

cat twelve_spp_proteins/longest_*.fa | seqkit rmdup - > makeblastdb_input.fa

第二步：All-by-all BLAST

makeblastdb -in makeblastdb_input.fa-dbtype prot -out blastdb

blastp -num_threads 30 -db blastdb -query makeblastdb_input.fa -outfmt 7 -seg yes > blast_output.txt

注：-seg参数过滤低复杂度的序列(即氨基酸编码为X),-num_threads线程数，此处设置为30。

第三步：使用MCL进行序列聚类

根据BLAST输出中序列相似性信息寻找聚类。这些聚类将是后续用于CAFE分析的基因家族。聚类这一步将通过mcl处理。使用shell命令将BLAST转成MCL能够识别的ABC格式（其实就是挑选三列，两个ID和Evalue）。

grep -v "#"  blast_output.txt | cut -f 1,2,11 >blast_output.abc

然后，创建网络文件(.mci)和字典文件(.tab)，然后进行聚类

mcxload -abc blast_output.abc --stream-mirror --stream-neg-log10 -stream-tf  'ceil(200)' -o blast_output.mci -write-tab blast_output.tab

其中：--stream-mirror: 为输入创建镜像，即每一个X-Y都有一个Y-X

          --stream-neg-log10: 对输入的数值做-log10转换

          -stream-tf: 对输入的数值进行一元函数转换，这里用到的是ceil(200)

根据mci文件进行聚类, 其中主要调整的参数是-I, 它决定了聚类的粒度，值越小那么聚类密度越大，这个值没有想象中的那么至关重要。一般设置为3，你也可以尝试用其他值，然后比较结果。最终的目的是正确分析物种间的直系同源基因。

mcl blast_output.mci -I 3

mcxdump -icl out.blast_output.mci.I30 -tabr blast_output.tab -o dump.blast_output.mci.I30

第四步：整理MCL的输出结果

上一步MCL的输出还不能直接用于CAFE，还需要对其进行解析并过滤。

 

第一步是将原来的mcl格式转成CAFE能用的格式。

pythonpython_scripts/cafetutorial_mcl2rawcafe.py  -i dump.blast_output.mci.I30 -o unfiltered_cafe_input.txt -sp "ENSG00 ENSPTR ENSPPY ENSPAN ENSNLE ENSMMUENSCJA ENSRNO ENSMUS ENSFCA ENSECA ENSBTA"

这里的"ENSG00" 是ENSEMBL编号中物种的标识符。

unfiltered_cafe_input.txt文件如下所示：

第二步，将那些基因拷贝数变异特别大的基因家族剔除掉，因为它会造成参数预测出错。下面的脚本是过滤掉一个或多个物种有超过100个基因拷贝的基因家族，虽然不是特别的严格，但效果和根据拷贝数变异过滤类似。

pythonpython_scripts/cafetutorial_clade_and_size_filter.py -iunfiltered_cafe_input.txt -o filtered_cafe_input.txt –s

然后把ID换成物种名字：

sed -i -e 's/ENSPAN/baboon/' -e 's/ENSFCA/cat/' -e 's/ENSBTA/cow/' -e's/ENSNLE/gibbon/' -e 's/ENSECA/horse/' -e 's/ENSG00/human/' -e's/ENSMMU/macaque/' -e 's/ENSCJA/marmoset/' -e 's/ENSMUS/mouse/' -e's/ENSPPY/orang/' -e 's/ENSRNO/rat/' -e 's/ENSPTR/chimp/' filtered_cafe_input.txt

sed -i -e 's/ENSPAN/baboon/' -e 's/ENSFCA/cat/' -e 's/ENSBTA/cow/' -e's/ENSNLE/gibbon/' -e 's/ENSECA/horse/' -e 's/ENSG00/human/' -e's/ENSMMU/macaque/' -e 's/ENSCJA/marmoset/' -e 's/ENSMUS/mouse/' -e's/ENSPPY/orang/' -e 's/ENSRNO/rat/' -e 's/ENSPTR/chimp/' large_filtered_cafe_input.txt

第五步：物种树推断

构建物种树主要分为多序列联配和系统发育树推测两步，之后在已有进化树的基础上构建超度量树用作CAFE输入。

 

多序列联配一般用的是单拷贝的直系同源基因（其实前面的OrthoFinder就生成的有），分别进行多序列联配之后然后合并成单个文件。接着用系统发育树推测软件进行建树，可选软件有

极大似然法: RAxML, PhyML, FastTree

贝叶斯法: MrBayes

推断超度量树

超度量树(ultrametric tree)也叫时间树，就是将系统发育树的标度改成时间，从根到所有物种的距离都相同。构建方法有很多，比较常用的就是r8s.

 

这里用cafetutorial_prep_r8s.py构建r8s的批量运行脚本，然后提取超度量树。

pythonpython_scripts/cafetutorial_prep_r8s.py -i twelve_spp_raxml_cat_tree_midpoint_rooted.txt -o r8s_ctl_file.txt -s 35157236 -p 'human,cat' -c '94'

/gpfs03/home/jingjing/software/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt

tail -n 1 r8s_tmp.txt | cut -c 16- >twelve_spp_r8s_ultrametric.txt

运行CAFE

运行CAFE有两种模式，一种是CAFE的命令行模式，先执行cafe进行CAFE的shell, 然后在其中执行命令。另一种是脚本模式，也就是你先把命令编辑完成，然后用cafe script_to_run.sh运行。

CAFE的主要功能就是根据给定的进化树和基因家族数估计一个或多个 birth-death()参数。参数描述的是基因出现或者消失的概率。

 

编辑cafetutorial_run1.sh。CAFE的命令不能有额外的空格出现在 tree后面的()中，以及lambda 的 -t 后的()中，否则运行时会无法正确解析文件导致报错。

mkdir -p reports

cafe cafetutorial_run1.sh

这步运行结束后的报告文件在reports/reportrun1.cafe，可以用已有的脚本分析哪些基因家族发生了扩张或者搜索。

python  /gpfs03/home/jingjing/software/CAFE/script/Fulton_python_scripts/cafetutorial_report_analysis.py  -i reports/report_run1.cafe -o reports/summary_run1   （注意这些python程序要基于python2才行）

在reports文件夹下会出现如下文件

 

summary_run1_node.txt: 统计每个节点中扩张，收缩的基因家族数

summary_run1_fams.txt: 具体发生变化的基因家族

看下CAFE的输出结果：

Lambda是整个进化树的预测值

 

# IDs of nodes表示不同节点的编号，这里cat为0，horse为2，cat和horse所在的节点是1.

 

最后是每个基因家族的结果。以最开始的表示行为例，第一列对应输入基因家族的编号；第二列是Newick的进化树，cat_61中的61表示该基因家族在cat里有61个基因；第三列是Family-wide P-value，用于表明该基因家族是否是显著性的扩张或是收缩，这里是0.124，说明变化不明显。在第三列的p值小于0.01时，第四列表明哪个分支的基因家族发生了变化，上图中只有ID 13的基因家族有变化。

cafe结果可视化

在网上搜cafe可视化发现并没有什么资料，都是说自己进行结果提取和画图。

这里有两种方法，一是使用cafe自带脚本计算出扩张和收缩数目后自行绘图，比如ggtree或者itol手动绘制等。

方法如下：

python cafetutorial_draw_tree.py -isummary_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，也比较麻烦。而且我自己测试时全部使用最新版并没有出现问题，所以安装时不需要降级。

 

但是运行的时候又报错了：

ERROR:qt.qpa.screen: QXcbConnection: Could not connect to display :0.0解决办法：

echo "exportQT_QPA_PLATFORM='offscreen'" >>~/.bashrc

source ~/.bashrc

python CAFE_fig.py example_result.cafe -pb 0.05 -pf 0.05 --dump test/ -g pdf --count_all_expansions  //CAFE_fig自带的例子

注意：CAFE_fig.py中对pixels_per_mya进行修正。源码中是1.0，如果想让长度变成5倍则更正为5.0即可。

其中一个家族的结果：

 

跑上面的例子的结果（然后由于未知原因，我画出的图片没有颜色，暂时还不清楚为什么）：

本文使用 文章同步助手 同步