基因家族扩张或收缩，进化树1 2

生物学知识：基因家族的形成方式 转载自https://www.omicsclass.com/article/494

基因家族（英语：Gene family），是来源于同一个祖先，由一个基因通过基因重复而产生两个或更多的拷贝而构成的一组基因，它们在结构和功能上具有明显的相似性，编码相似的蛋白质产物。

基因家族的形成（加倍和扩张）

1.染色体組多倍化( Whole genome duplication or polyploidization)（WGD）：大部分植物都经历过古老的全基因组复制事件（WGD）or多倍化，它是一种大规模的染色体倍增事件，一次性增加一个物种所有基因的剂量，导致基因组中保留着大量染色体倍增的片段。全基因组复制造成的重复区通常是一大片区域中所有基因的重复，而不是单个基因或几个基因的重复。高等植物，很多都经历过多倍化过程如棉花（AADD），花生（AABB），小麦（AABBDD）等等，那么基因组在多倍化的过程中，就会发生基因的成倍增加：
2.串联复制(Tandem duplication) (TD)：串联复制主要发生在染色体重组区域，串联复制形成的基因家族成员通常紧密排列在同一条染色体上，形成一个序列相似、功能相近的基因簇。
3.染色体片段复制(Segmental duplication)(SD)：导致复制的基因距离较远，甚至位于不同的染色体。
4.反转录转座(Retrotransposition)：反转录转座是指件已经转录和剪切的mRNA，再经过逆转录过程形成cDNA，然后随机插入到染色体的某一位置形成新的重复基因的过程。经反转录转座形成的新基因往往由于缺少必要的调控序列，通常都是不能表达的假基因。
5.外显子重复和改组（Exon duplication and shuffling ）
大片段复制的家族基因，及反转录转座也会导致加倍的家族基因位于不同的染色体
严格区分的话，应该用共线性的方法去查找大片段复制的基因，如mcscan去查找：因为大片段复制的基因周围的基因也有相似性，而反转录转座就不同了，周围的其他基因不相似。

0.收集需要比较的物种的蛋白序列，并获取每个基因编码的最长的蛋白序列

1.使用OrthoFinder查找直系同源基因

2.Iqtree进化树构建

3.MCMCTREE估计分化时间，R8S用推断的物种树提取超度量树（ultrametric tree）

R8s的下载安装

wget -q https://sourceforge.net/projects/r8s/files/r8s1.81.tar.gz /
  && tar -zxvf r8s1.81.tar.gz /
  && cd r8s1.81/src /
  && cp Makefile.linux Makefile.linux.bak /
  && sed -i 's/continuousML.o //' Makefile.linux /
  && sed -i 's/continuousML.o:/#continuousML.o:/' Makefile.linux /
  && make -f Makefile.linux

4. 基因家族扩张与收缩分析（CAFE）

cafe结果可视化

安装cafe_fig，用于cafe产出文件可视化。

pip3 install ete3
git clone https://github.com/LKremer/CAFE_fig.git
cd CAFE_fig
python3 CAFE_fig.py -h

安装如果报错，缺少依赖包，就用pip3安装。
如果报错NameError: name 'ListType' is not defined，pip3 install scipy可解决。
错误：qt.qpa.plugin: Could not load the Qt platform plugin "xcb" in "" even though it was found
解决办法：用conda安装pyqt5 conda install pyqt
ERROR:qt.qpa.screen: QXcbConnection: Could not connect to display :0.0
解决办法：

echo "export QT_QPA_PLATFORM='offscreen'" >>~/.bashrc
source ~/.bashrc

运行例子

python3 CAFE_fig.py example_result.cafe -pb 0.05 -pf 0.05 --dump test/ -g svg --count_all_expansions
示例输出结果，需要自己用AI调整SVG.

image.png

cafe结果可视化方法2：利用caft自带的脚本可视化

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

cafe自带的py脚本在其他人的github

cafe分析的准备

awk -F"/t" '{print "Desc""/t"$1"/t"$2"/t"$3"/t"$4"/t"$5"/t"$6"/t"$7"/t"$8"/t"$9"/t"$10"/t"$11"/t"$12"/t"$13}' Orthogroups.GeneCount.tsv > cafe-input.data 
sed -e '1 s/(null)/Desc/' -e '1 s/Orthogroup/Family/' -i cafe-input.data
python cafetutorial_clade_and_size_filter.py -i cafe-input.data -o filtered_cafe_input.txt -s 2> filtered.log 

生成cafe脚本程序

#!/share/home/chaimao1/miniconda3/bin/cafe
version
date

load -i ortho_input.txt -t 8 -p 0.01
tree 

lambda -s
report out' >cafe_command
chmod 757 cafe_command

caferror.py -i cafe_command

cafe程序的命令参数简介

version
显示版本信息

date
显示当前日期

load
-i 输入的数据文件ortho_input.txt ,是基因家族在各物种中的数目，该文件有多列，以tab分割。第一行格式固定
-t 运行的线程数量，默认是8
-l 设置程序输出的日志文件
-p 设置p-value的阈值，默认是0.01

tree (((chinmp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93) 输入超度量树
lambda
-l 设置λ的值
-s 自动寻找最优的λ值
-t 默认情况下所有的分支的λ值是相同的，如果需要设置不同的λ值。则该参数应该和上面的树文件格式一样。(((1,1)1,(2,2)2)2,2)

report out
设置输出文件的前缀为out

ortho_input.txt是统计所有物种的基因家族的基因的数量，可以从Orthofinder的输出结果Results_Jul08/Orthogroups/Orthogroups.GeneCount.tsv中获取。
ortho_input.txt的文件头部格式必须是下面这种，中间的分割符是tab.

Description  ID  物种1  物种2   ...  物种n

awk '{print Description"/t"$0 }' Orthogroups.GeneCount.tsv|sed -e '1 s/Orthogroup/Family/g' >cafe-input.data
python cafetutorial_clade_and_size_filter.py -i cafe-input.data -o filtered_cafe_input.txt -s 2> filtered.log