WeFun blog
0%

基因家族扩张收缩

发表于
本文字数： 10k 阅读时长 ≈ 9 分钟

基因家族指来自大部份物种的MRCA（最近共同祖先，Most Recent Common Ancestor）的同一个始祖基因演化而来的一组基因。
CAFE(Computational Analysis of gene Family Evolution)是一款以解释系统发育历史的方式分析基因家族大小变化的软件，这种分析常被称为基因家族收缩扩张(Gene family expansions and contractions)分析。2005年该软件相关算法就已经发表，最近一次更新在2023年。CAFE使用出生和死亡过程来模拟用户指定的系统发育树中的基因获得和丢失（获得替换率，λ），2017年时首次提出不同枝基因的出生死亡率不同，这样可以算出父节点到子节点的基因家族大小转移率，也可推断祖先物种的基因家族大小

研究基因家族在进化过程中的大小变化，即关注旁系同源基因的演化，本次教程主要来源： https://yanzhongsino.github.io/2021/10/29/bioinfo_gene.family_CAFE5/

1.安装

1
conda install -c bioconda cafe

安装的是CAFE4
直接上官网安装可以安装到最新版

1
2
3
4
git clone https://github.com/hahnlab/CAFE5.git
cd CAFE5
./configure
make

如果遇到报错

1
2
3
4
5
6
g++ -std=c++11 -I. -O3 -include config.h -DDOCTEST_CONFIG_DISABLE -fopenmp  -c src/matrix_cache.cpp -o obj/matrix_cache.o
src/matrix_cache.cpp:12:10: fatal error: cblas.h: No such file or directory
   12 | #include "cblas.h"
      |          ^~~~~~~~~
compilation terminated.
make: *** [Makefile:35: obj/matrix_cache.o] Error 1

这是因为缺少数学运算程序，只需要安装blas和lapack即可，主要教程为： https://zhuanlan.zhihu.com/p/520848641
blas已经包括在lapack中，因此只需要下载lapack即可 

1
2
#下载LAPACK
wget https://github.com/Reference-LAPACK/lapack/archive/refs/tags/v3.12.1.tar.gz

解压完成之后进入lapack-3.12.1

1
2
3
4
5
6
cd lapack-3.12.1
cp make.inc.example make.inc
make blaslib
make cblaslib
make lapacklib
make lapackelib

编译完成之后将静态库加入到环境变量即可（没有root权限）

1
2
export LIBRARY_PATH="$LIBRARY_PATH:/xxx/lapack-3.10.1"
export C_INCLUDE_PATH="$C_INCLUDE_PATH:/xxx/lapack-3.10.1/LAPACKE/include:/xxx/lapack-3.10.1/CBLAS/include"

除此之外，还需要将CAFE5/makefile文件开头部分加上cblas.h位置

1
CFLAGS += -I/home/walnut168/mnt/wf/program/lapack-3.12.1/CBLAS/include

2.准备文件

CAFE5需要至少两个输入文件，一个是基因家族计数文件gene_families.txt，一个是树文件tree.txt

2.1gene_families.txt

制表符分隔的基因家族计数文件，通常用OrthoMCL, SwiftOrtho, FastOrtho, OrthAgogue, OrthoFinder等软件获取，因此可以直接使用OrthoFinder的Orthogroups.GeneCount.tsv文件来获得
文件位置为/data/OrthoFinder/Results_Apr21/Orthogroups/Orthogroups.GeneCount.tsv

CAFE给的示例文件文件如下：

1
2
3
4
5
6
7
#第一列是功能描述Desc，不清楚可以写为null，第二列是Family ID
Desc	Family ID	human	chimp	orang	baboon	gibbon	macaque	marmoset rat	mouse	cat	horse	cow
ATPase	ORTHOMCL1	 52	 55	 54	 57	 54	  56	  56	 53	 52	57	55	 54
(null)	ORTHOMCL2	 76	 51	 41	 39	 45	  36	  37	 67	 79	37	41	 49
HMG box	ORTHOMCL3	 50	 49	 48	 48	 46	  49	  48	 55	 52	51	47	 55
(null)	ORTHOMCL4	 43	 43	 47	 53	 44	  47	  46	 59	 58	51	50	 55
Dynamin	ORTHOMCL5	 43	 40	 43	 44	 31	  46	  33	 79	 70	43	49	 50

使用Orthogroups.GeneCount.tsv得到gene_families.txt，需要进行一些处理

1
awk -v OFS="\t" '{$NF=null;print $1,$0}' Orthogroups.GeneCount.tsv |sed -E -e 's/Orthogroup/desc/' -e 's/_[^\t]+//g' >gene_families.txt

这一步的主要操作是给文件加上第一列desc，除此之外就是将文件分隔符改为\t

值得注意的是如果不事先删除不同物种间的基因拷贝数变异特别大的基因家族，则会出现报错 

Failed to initialize any reasonable values

可以使用cafe5自带的一个脚本，cafe5的tutorial目录下脚本clade_and_size_filter.py可以筛除一个或以上物种有超过100个基因拷贝的基因家族。

1
python clade_and_size_filter.py -s -i gene_families.txt -o gene_familie_filter.txt

注意运行脚本时gene_families.txt需删去表头，但是CAFE的输入文件需要表头，因此后面需要加上

2.2tree.txt

带分化时间的，二叉的，有根的，超度量树，树的格式为newick格式，超度量树指树的最近发散序列的分支长度相等。
可以使用PAML生成的结果FigTree.tre转换成tree.txt文件 

1
2
3
grep -E -v "NEXUS|BEGIN|END" FigTree.tre|sed -E -e "s/\[[^]]*\]//g" -e "s/[ \t]//g" -e "/^$/d" -e "s/UTREE/tree tree/" >tree.txt
cat FigTree.tre | sed 's/\[[^]]\+\]//g'| awk -F "=" '/UTREE/{print $2} ' > input.tree.nwk
sed -e 's/:/\n:/g' -e 's/\([),]\)/\n\1/g' input.tree.nwk | awk '{if($1~/:$/){printf ":"100*$2} else {printf $0}}' | sed 's/\s\+//g' > input.tree

3.运行

1
2
3
4
5
6
cafe5 --infile gene_families.txt --tree tree.txt --cores 5 --output_prefix cafe_out 

--infile gene_families.txt 输入基因家族文件
--tree tree.txt 输入树文件
--cores 5 指定5个核运行
--output_prefix cafe_out 输出文件名

其他可添加参数

  • -p代表指定root frequency distribution为泊松分布（默认是均匀分布uniform distribution）。
  • -k 3代表使用GAMMA模型（默认是base模型）并且使用3种gamma rate（代表不同基因家族有着不同的进化速率）。-k的值需要运行多次比较likelihood并确保收敛后才知道使用哪个最好，一般来说2-5之间试一试。

4.结果

1
2
3
4
5
6
7
8
9
cafe_out/
Base_asr.tre # 进化树格式，记录了基因数目, 显著变化节点*标
Base_branch_probabilities.tab # 记录每个节点pvalue
Base_change.tab # 记录每个节点增加或减少的基因数目
Base_clade_results.txt # 每个节点收缩扩张基因家族数目统计，可用于做饼图
Base_count.tab # 记录每个节点基因数目
Base_family_likelihoods.txt
Base_family_results.txt # 记录家族水平是否显著变化即p<0.05
Base_results.txt

Base_asr.tre文件是记录信息最全的文件，开头部分如下：

1
2
3
4
#nexus
BEGIN TREES;
  TREE OG0000000 = ((LZ<1>*_3:8.21782,WX<2>_34:8.21782)<27>*_43:1.706,(JB<3>_69:9.62703,(SJ<4>_72:8.85394,(YX<5>*_116:7.83667,(WYZ<6>*_49:6.95006,(SL<7>_65:6.43522,(SD<8>_88:5.94193,((DX<9>*_188:4.89218,(PT<10>_112:4.47741,((MGY<11>*_144:2.19967,NNJ<12>*_102:2.19967)<38>_117:1.19645,(YJ<13>*_96:2.48949,(nigra<14>_134:0.284379,CH<15>_137:0.284379)<40>*_134:2.20511)<39>_113:0.906634)<37>_112:1.08128)<36>_107:0.414774)<35>*_105:0.614199,(MDL<16>*_58:4.33793,((EL<17>_95:1.21337,XZ<18>_90:1.21337)<43>_90:1.9202,(EZQ<19>*_62:1.90236,(HG<20>_76:1.32269,((TN<21>*_85:1.16139,HYM<22>*_66:1.16139)<47>_77:0.072474,(XH<23>_82:0.435824,(WC<24>_81:0.35856,JN<25>_80:0.35856)<49>_80:0.077264)<48>_80:0.798037)<46>_77:0.088832)<45>_77:0.579668)<44>_76:1.23121)<42>_82:1.20436)<41>*_83:1.16845)<34>*_91:0.43555)<33>*_86:0.493286)<32>_81:0.514848)<31>_77:0.886606)<30>*_74:1.01727)<29>*_66:0.773085)<28>_60:0.296801)<26>_58;
  TREE OG0000005 = ((LZ<1>*_1:8.21782,WX<2>_12:8.21782)<27>*_18:1.706,(JB<3>*_41:9.62703,(SJ<4>*_13:8.85394,(YX<5>_32:7.83667,(WYZ<6>_49:6.95006,(SL<7>*_32:6.43522,(SD<8>_60:5.94193,((DX<9>*_72:4.89218,(PT<10>*_84:4.47741,((MGY<11>_63:2.19967,NNJ<12>*_45:2.19967)<38>_55:1.19645,(YJ<13>_60:2.48949,(nigra<14>_42:0.284379,CH<15>_44:0.284379)<40>*_44:2.20511)<39>_54:0.906634)<37>_56:1.08128)<36>_59:0.414774)<35>_58:0.614199,(MDL<16>*_45:4.33793,((EL<17>_73:1.21337,XZ<18>*_85:1.21337)<43>_77:1.9202,(EZQ<19>*_94:1.90236,(HG<20>*_88:1.32269,((TN<21>_70:1.16139,HYM<22>_79:1.16139)<47>_74:0.072474,(XH<23>_64:0.435824,(WC<24>*_48:0.35856,JN<25>_58:0.35856)<49>*_58:0.077264)<48>*_61:0.798037)<46>_74:0.088832)<45>_75:0.579668)<44>_77:1.23121)<42>*_71:1.20436)<41>_60:1.16845)<34>*_55:0.43555)<33>*_51:0.493286)<32>_46:0.514848)<31>*_43:0.886606)<30>*_36:1.01727)<29>_29:0.773085)<28>_27:0.296801)<26>_26;

Base_asr.tre记录文件信息为TREE <Family ID> +树(物种名<内部节点名称>_基因家族数量：枝长……..以此类推)，*号表示该节点相对于最近祖先的基因家族变化非常显著，p<0.05

5.可视化

5.1绘制饼图

将基因家族扩张收缩的条数用于绘制饼图，绘制过程可以使用itol网站，绘制文件是树文件，因此需要从CAFE结果中提取树文件
Base_asr.tre中记录了所有基因家族的树，所有的树都是一致的，只不过基因家族数不同，因此可以使用以下代码提取任意一个树即可 

1
awk '/TREE OG0000000 = /{print $NF}' out/Base_asr.tre |sed -r 's/_[0-9]+//g ; s/([^)])<[0-9]+>/\1/g ; s/\*//g' >out.tree

然后就是具体基因家族扩张收缩绘图文件获得

1
2
3
4
5
6
7
8
9
10
11
##提取cafe输出的家族总个数，用于计算基因家族有多少基因家族没有变化
Tnum=$(grep -c "OG" Orthogroups.GeneCount.tsv)
##生成用于itol绘图的dataset文件
echo -e "DATASET_PIECHART
SEPARATOR\tTAB
DATASET_LABEL\ttest
FIELD_COLORS\t#339900\t#FF0033\t#0099FF
FIELD_LABELS\tIncrease\tDecrease\tRemain
DATA" > cafe_itol.txt
##获得扩张收缩家族数
sed -r 's/^([^>]+)<[0-9]+>/\1/g' out/Base_clade_results.txt | awk  '!/^#/{ p=0.5; s=10; C='$Tnum'-$2-$3-1; if(!/^</){p=-1}; print $1 "\t" p "\t" s "\t" $2 "\t" $3 "\t" C }' >> cafe_itol.txt

itol的dataset文件表头以下信息的含义为
第一列绘制饼图的内部节点的名称，第二列为饼图绘制的位置，第四列为饼图的直径，第五列之后是绘制饼图的绘制数据
绘制时，直接拖入饼图dataset文件即可

除此之外，itol网站还可以添加标签注释，与上面类似，只需要写好绘图的dataset文件即可，所有文件格式均在itol最上面HELP->HELP PAGE，页面左侧Dataset type中有详细介绍与示例文件，以下展示添加文本标签的示例文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
DATASET_TEXT
# 在文本数据集中，每个ID关联一个文本标签，可直接显示在节点分支上或树外
# 以井号(#)开头的行是注释，解析时会被忽略
#=================================================================#
#                        必填设置                                 #
#=================================================================#
# 选择下方数据的分隔符（TAB, SPACE 或 COMMA）。此分隔符必须在全文件统一使用。
# SEPARATOR TAB
# SEPARATOR SPACE
SEPARATOR COMMA  # 当前使用逗号分隔

# 用于图例表的标签（后续可修改）
DATASET_LABEL,example text dataset

# 数据集颜色（后续可修改）
COLOR,#ff0000  # 红色

#=================================================================#
#                        可选设置                                 #
#=================================================================#
# 所有其他可选设置可在网页界面中修改（位于'Datasets'选项卡下）

# 左边距（可正可负），用于调整与下一个数据集的间距。仅影响外部文本标签
MARGIN,0

# 仅适用于外部文本标签。设为1时，即使内部节点未折叠也会显示其标签（可能导致重叠）
SHOW_INTERNAL,0

# 所有标签旋转指定角度
ALL_LABELS_ROTATION,0

# 默认内部标签在分支上方。设为1时，标签显示在分支下方
LABELS_BELOW,1

# 垂直方向移动标签的像素值（可正可负）
VERTICAL_SHIFT,0

# 设为1时，树形旋转不会影响单个标签的旋转
STRAIGHT_LABELS,0

# 仅适用于外部标签。设为1时，标签将沿树对齐（圆形模式）或垂直显示（常规模式），忽略所有旋转参数
ALIGN_TO_TREE,0

# 字体大小因子：外部标签默认字体大小略小于叶节点间距，此值可缩放（<1减小，>1增大）
SIZE_FACTOR,1

# 外部标签额外水平偏移（适用于无根模式调整标签位置）
EXTERNAL_LABEL_SHIFT,0


# 是否在外部标签列上方显示数据集标签（注释状态）
#SHOW_LABELS 1

# 数据集标签大小因子（注释状态）
#LABEL_SIZE_FACTOR 1

# 数据集标签旋转（注释状态）
#LABEL_ROTATION 0

# 数据集标签像素偏移（注释状态）
#LABEL_SHIFT 0

# 数据集标签与树对齐（仅圆形模式有效，注释状态）
#LABEL_ALIGN_TO_TREE,0

# 内部节点可通过ID直接指定，或用iTOL帮助文档描述的'最后共同祖先'方法
#=================================================================#
#                 数据部分（"DATA"关键字后开始）                  #
#=================================================================#
# 每个节点可包含以下字段：
# ID,标签文本,位置,颜色,样式,大小因子,旋转角度

# 位置参数说明：
#  -1 = 外部标签
#  0-1 = 内部标签（0:分支起点, 0.5:中点, 1:终点）
# 样式可选：'normal','bold','italic','bold-italic'
# 大小因子将乘以标准字体大小

DATA
# 示例：

# 节点9598将显示红色粗体外部标签"Pan troglodytes"，字体大小为标准2倍
#9598,Pan troglodytes,-1,#ff0000,bold,2,0

# 节点9606将显示多样式混合外部标签
#9606,<bi color='#006600'>Homo </bi><i>sapiens</i><sup size='0.5' color='#999999'>citation</sup>,-1,#000000,normal,1,0

# 节点4530将显示蓝色粗斜体内部标签"Oryza sativa"，位于分支起点
#4530,Oryza sativa,0,#0000ff,bold-italic,1,0

其实与上面饼图文件一样，表头信息加后面的数据
echo -e "DATASET_TEXT
SEPARATOR	TAB
DATASET_LABEL	TEST
COLOR	#ff0000
DATA" > cafe.txt
sed -r 's/^([^>]+)<[0-9]+>/\1/g' |awk '{print $1"\t+"$2"/-"$3"\t0.5\t#ff0000\tnormal\t1\t0"}' >> cafe.txt
sed -i '/<[0-9]\+>/d' cafe.txt

配置完成之后一样，拖入网站即可