教程 | 基因家族扩张与收缩分析

介绍

#1

■ 在比较基因组分析，对直系同源基因分析后，往往接着就是基因家族扩张收缩分析。确定生物间表型差异背后的遗传变化和导致变化的进化压力，是进化生物学的主要目标之一。基因组分析工作已经揭示了物种间基因家族的成员频繁获得和丢失。基因家族规模的变化可能有利、有害或者中性，但基因家族的数量变化也是形成物种特异的重要原因之一。

#2

■ 推荐一下普遍使用的软件CAFE (Computational Analysis of gene Family Evolution)。小编永远支持具体问题还是具体分析，生物分析软件也只是辅助我们获得一个预测性的结果。

#3

■首先，CAFE由Matthew W. Hahn 课题组在2005年提出评估基因家族进化速度和模式的模型，2006年CAFE软件面世，2013年推出CAFE 3软件，2020年更新CAFE 5软件。（可见该软件还是一直更新）

#4

■ 摘一段CAFE5（Bie et. al , 2021）简短解释CAFE工作原理：

The probabilistic model adopted in CAFE was introduced by Hahn et al. (2005); it uses a random birth and death process to model gene gain and loss along each lineage of a phylogenetic tree. In order to make inferences over a whole phylogeny, a probabilistic graphical model (Lauritzen, 1996; M. I. Jordan, manuscript in preparation) is used to calculate the probability of transitions in gene family size from parent to child nodes in the phylogeny.Using the graphical models machinery, one can draw inferences on the gene family size for all ancestral species.

#5

■ CAFE 应用随机出生死亡的模型，模拟一个系统发育过程基因家族得失。为了推断系统发育过程，可计算由父节点到子节点的基因家族大小转移率，也可推断祖先物种的基因家族大小。

安装

#CAFE5下载
git clone https://github.com/hahnlab/CAFE5.git
cd CAFE5/
mkdir -p `pwd`/bin
make install prefix=`pwd`

准备文件

需要Orthofinder的结果文件

可以再翻查以前的推文
Orthofinder下篇
Orthofinder上篇

mkdir CAFE && cd CAFE
#定义我Orthofinder结果路径
OrthofinderPATH=$mypath/OrthoFinder/Results_Jun12
cp $OrthofinderPATH/Orthogroups/Orthogroups.GeneCount.tsv .

需要时间分歧树文件

也可以翻查过去的mcmctree推文

cp $mcmctreePATH/FigTree.tre .

整理输入文件

GeneCounts.tsv

#查看Orthofinder结果文件
head Orthogroups.GeneCount.tsv
#Orthogroup      Cpapaya Csinensis       Dzibethinus     Fvesca  Tcacao  Vvinifera       Total
#OG0000000       374     6       0       3       1       0       384
#OG0000001       308     0       1       1       15      6       331
#并不符合CAFE输入的文件格式，CAFE似乎是默认格式和OrthoMCL一致，所以需要调整
#
awk -v  OFS="\t" '{if($1=="Orthogroup"){print"Descript",$1,$2,$3,$4,$5,$6,$7}else{print"Null",$1,$2,$3,$4,$5,$6,$7}}' Orthogroups.GeneCount.tsv > GeneCounts.tsv

查看分歧树
需要分歧树的节点时间，并不需要置信区间。

#NEXUS
BEGIN TREES;

        UTREE 1 = (Vvi: 1.017597, (Fve: 1.012650, (Csi: 0.888532, ((Dzi: 0.267560, Tca: 0.267560) [&95%={0.241975, 0.297312}]: 0.517603, Cpa: 0.785163) [&95%={0.736974, 0.849596}]: 0.103369) [&95%={0.839088, 0.957735}]: 0.124118) [&95%={0.963919, 1.08662}]: 0.004947) [&95%={0.968735, 1.09186}];

END;

runcafe.bash
输入的树文本内不要又空格

#!cafe
load -i GeneCounts.tsv  -t 20 -l log.txt
tree (Vvinifera:1.017597,(Fvesca:1.012650,(Csinensis:0.888532,((Dzibethinus:0.267560,Tcacao:0.267560):0.517603,Cpapaya:0.785163):0.103369):0.124118):0.004947)ivaana;
lambda -s -t (1,(1,(1,((1,1)1,1)1)1)1)
report out

load参数
-i 输入的数据文件
-t 设置程序运行的线程数，默认为 8
-l 设置输出的日志文件，默认标准输出
-p 设置 p_value 的阈值，默认为 0.01

运行

cafe runcafe.bash
#生成文件
#out.cafe
head out.cafe

查看文件内容

#第一列是使用的时间分歧树内容
Tree:(Vvinifera:1.0176,(Fvesca:1.01265,(Csinensis:0.888532,((Dzibethinus:0.26756,Tcacao:0.26756):0.517603,Cpapaya:0.785163):0.103369):0.124118):0.004947)ivaana
#推断的出生-死亡参数λ预测值
Lambda: 0.384269
Lambda tree:    (1,(1,(1,((1,1)1,1)1)1)1)
#每个节点编号
# IDs of nodes:(Vvinifera<0>,(Fvesca<2>,(Csinensis<4>,((Dzibethinus<6>,Tcacao<8>)<7>,Cpapaya<10>)<9>)<5>)<3>)ivaana<1>
# Output format for: ' Average Expansion', 'Expansions', 'No Change', 'Contractions', and 'Branch-specific P-values' = (node ID, node ID): (0,3) (2,5) (4,9) (6,8) (7,10)
# Output format for 'Branch cutting P-values' and 'Likelihood Ratio Test': (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
#平均每个基因家族中扩张的基因数目，负数表示基因家族收缩
Average Expansion:      (-0.0456955,-0.0292223) (1.53621,0)     (0.201884,0)    (0.446966,0.0637459)    (0,0.000175246)
#发生了扩张的基因家族数目
Expansion :     (3777,0)        (11643,0)       (3901,0)        (8934,3795)     (0,3380)
#没有发生改变的基因家族数目
nRemain :       (10365,22158)   (6051,22825)    (11427,22825)   (6762,12401)    (22825,11162)
#发生了收缩的基因家族数目
nDecrease :     (8683,667)      (5131,0)        (7497,0)        (7129,6629)     (0,8283)
#每个基因家族具体扩张和收缩的情况，ID、树的信息（物种下划线后是该基因家族的数量）、基因家族扩张收缩的pvalue、每个分枝的基因家族扩张的显著性
'ID'    'Newick'        'Family-wide P-value'   'Viterbi P-values'      'cut P-value'   'Likelihood Ratio'
OG0000000       (Vvinifera_0:1.0176,(Fvesca_3:1.01265,(Csinensis_6:0.888532,((Dzibethinus_0:0.26756,Tcacao_1:0.26756)_0:0.517603,Cpapaya_374:0.785163)_0:0.103369)_0:0.124118)_0:0.004947)ivaana_1      0.5005  ((-,-),(-,-),(-,-),(-,-),(-,-))

提取更详细的扩张收缩信息
使用安装包里的python脚本。如果是conda安装的，可以从git-hub中下载该部分内容；

#使用python2运行，如果环境没python2可以安装。
python2 $CafePATH/python_scripts/cafetutorial_report_analysis.py -i out.cafe  -o out.report
#形成文件
#out.report_anc.txt #和GeneCounts.tsv内容相似，格式不一样
#out.report_fams.txt #每个节点发生扩张收缩的基因家族，有*是具有显著性。
#out.report_pub.txt #每个物种收缩、扩张、获得、丢失的情况
#out.report_node.txt #每个节点基因家族收缩、扩张和快速进化的家族数

画树

#也是使用安装包内的python脚本
#但如果使用的不是具有界面化的服务器，建议按照以下操作修改脚本
vi $CafePATH/python_scripts/cafetutorial_draw_tree.py
## 在import 模块后
## 添加一行不输出图片在桌面
plt.switch_backend('agg')

#运行，也是用python2
python2 $CafePATH/python_scripts/cafetutorial_draw_tree.py \
-i out.report_node.txt  -y Expansions  \
-t '(Vvinifera:1.0176,(Fvesca:1.01265,(Csinensis:0.888532,((Dzibethinus:0.26756,Tcacao:0.26756):0.517603,Cpapaya:0.785163):0.103369):0.124118):0.004947)'  \
-d '(Vvinifera<0>,(Fvesca<2>,(Csinensis<4>,((Dzibethinus<6>,Tcacao<8>)<7>,Cpapaya<10>)<9>)<5>)<3>)<1>' \
-o Expansionstree.png

-i 输入的信息文件
-y是对应输入文件的标题选择展示 可选：Expansions/Contractions/Rapid
-t 是输入树文件 （在cafe输出文件能找到）
-d 是树的结构文件（在cafe输出文件能找到）
-o 命名，脚本默认是png；
不过由于是由python来写，还是能修改脚本的命令调整输出为svg

#将这句的格式调整svg输出或者pdf都可
fig.savefig(output_file, format='svg', bbox_inches='tight', dpi=300)

输出的图展示

始终还是输出的不够好看（无法直接发文章），不过没关系总有方法进行美观的。