进化树的构建与美化

构建进化树常常有多种办法，常用的有邻接法（NJ），UMPGA，最大简约法（MP），最大似然法（ML），贝叶斯法，不同方法各有其优缺点，最常用的就是ML和贝叶斯法，各个方法的原理及优缺点不在此处进行赘述，以下仅介绍建树的具体流程
构建进化树一般分为两步，第一步进行多序列比对，第二步对比对结果建树。

1.多序列比对

多序列比对常用多序列比对常用软件为 muscle 和 mafft, 这里使用 muscle。github下载链接为https://github.com/rcedgar/muscle/archive/refs/tags/v5.3.tar.gz，官网位置为https://github.com/rcedgar/muscle

1 2	muscle -in all_FAD.fasta -out all_FAD.muscle.fasta #all_FAD.fasta为建树物种共同的某一基因家族的序列组合，为蛋白质序列

如果想要使用mafft可以使用以下命令：

#高速模式
mafft in > out
#高准确性模式（序列<200，aa/nt<2000）
mafft --maxiterate 1000 --localpair in > out
mafft --maxiterate 1000 --genafpair in > out
mafft --maxiterate 1000 --globarpair in > out

2.构建进化树

RAxML 是一款基于最大似然法构建系统发育树的软件，运算速度较快，推荐使用。RAxML 的输入文件为多序列比对文件，phylip 格式或者 fasta格式，输出文件为 newick 格式进化树
软件下载位置：https://cme.h-its.org/exelixis/web/software/raxml/hands_on.html
下载好后可以看到，有raxmlHPC和raxmlHPC-PHREADS 分别为单线程和多线程版本，一般使用后者

raxmlHPC-PTHREADS -m PROTGAMMAJTT -f a -p 123 -x 123 -# 100 -n out -s all_FAD.muscle.fasta 1>tree.log 2>tree.err

-m PROTGAMMAJTT \ # 指定建树模型，PROT表示输入为蛋白质，GAMMA表示离散位点使用的估计方式，JTT表示氨基酸模型
-f a \ # 指定算法，a表示使用rapidBootstrap方法并搜索最佳ML树
-p 123 \ # parsimonyRandomSeed，随机数种子
-x 123 \ # rapidBoostrap时RandomNumberSeed，随机数种子
-# 100 \ # bootstrap次数
-n out \ # 指定输出文件后缀
-s all_FAD.muscle.fasta \ # 指定输入的比对文件
1>tree.log 2>tree.err#存储日志
-T增加线程数，默认为2

#使用iqtree构建
iqtree 可以自动估计最佳模型，并进行 bootstrap 分析，也是一款推荐的建树软件。软件网址：http://www.iqtree.org/

iqtree -s all_FAD.muscle.fasta -m MFP -bb 1000 -nt 4-pre iqtree.all_FAD.muscle.fasta
-s all_FAD.muscle.fasta \ # 指定输入的比对文件
-m MFP \ # iqtree自动选择模型
-bb 1000 \ # bootstrap分析次数，至少是1000
-nt 4 \ # 线程数
-pre iqtree.all_FAD.muscle.fasta # 输出文件前缀

通常情况下我们无法直接获得建树物种共同的某一基因家族的序列的组合，此时我们可以借助之前orthofinder的结果来进行构建物种树，ortherfinder操作具体教程详见《查找共有和独有基因》

构建树的方式有两种，一种是使用连接法/串联法，将某一物种的所有单拷贝基因基因家族连接成一整条序列（supergene），每一物种都进行类似操作之后进行建树
另一种是合并法，每一个基因家族单独建树，总结其中出现次数最多的建树情况，然后将这种情况当作是最终结果

使用连接法时可以将pep文件转成CDS序列文件，再进行序列比对，进而进行建树，这种方式通常用于亲缘关系较近的物种建树

还可以使用4d位点，4d位点指四重简并位点，该一位点不论是ATCG中的任意碱基都不影响该蛋白质的编码，仅保留4d位点的序列可以减少比对时的数据量，进而提升建树速度及准确度

orthofinder使用msa的建树使用的就是连接法，不加msa使用的就是合并法，通常情况下建议使用msa参数

假设单拷贝基因家族的数量太少，此时可以将条件放宽，认为不是所有物种都具有的单拷贝基因家族才为共有，即假设共有10个物种，9个物种具有该基因家族，仍认为这是共有单拷贝基因家族

串联法
单拷贝基因家族位置为：data/OrthoFinder/Results_May14/MultipleSequenceAlignments/SpeciesTreeAlignment.fa

如果某一些位置的比对结果较差，可以选择将较差的部分删掉，可以使用trimal，基因树一般需要进行此操作，物种树一般不需要进行此操作

1
2
3

git clone https://github.com/scapella/trimal.git
cd trimal/source
make

trimal -in SpeciesTreeAlignment.fa -out SpeciesTreeAlignment_trim.fa -fasta -gt 0.6 -cons 60 

-gt 0.6 #保留60%基因共有的位点
-cons 60 # 过滤后总长不少于输入数据的60%

-fasta可替换
-CUSTAL格式的-clustal输出文件
-nbrf NBRF / PIR格式的输出文件
-nexus NEXUS格式的输出文件
-mega MEGA格式的输出文件
-phylip3.2 PHYLIP3.2格式的输出文件
-phylip PHYLIP / PHYLIP4格式的输出文件

1	raxmlHPC-PTHREADS -m PROTGAMMAJTT -f a -p 123 -x 123 -# 100 -n out -T 20 -s SpeciesTreeAlignment.fa 1>tree.log 2>tree.err

参数与上一致，不过多赘述

合并法
orthofinder已自建各基因家族树，位置为data/OrthoFinder/Results_Dec03/Orthogroups/Orthogroups_SingleCopyOrthologues.txt

先安装Astral

1	git clone https://github.com/smirarab/ASTRAL.git

#先使用脚本提取单拷贝基因树
orthDir=data3/OrthoFinder/Results_Dec03/
cat $orthDir/Orthogroups/Orthogroups_SingleCopyOrthologues.txt | \
while read aa; \
do cat $orthDir/Gene_Trees/$aa\_tree.txt | \
awk '{print $0}' ;\
done > SingleCopy.trees

#去除基因ID，保留物种名
sed 's/\([(,][A-Za-z]\+\)_[^:]*/\1/g' SingleCopy.trees > Astral_input.trees

#运行Astral
java -jar astral.5.7.8.jar -i Astral_input.trees -o Astral_output.tree 2>out.log

如果物种树较多，astral会出现内存不够的情况，实际测试时使用-Xmx300G添加了300G内存仍显示内存不足

如果Orthogroups_SingleCopyOrthologues.txt为空，需要放宽标准，选取90%或者95%的物种有的单拷贝同源序列可以使用脚本：

import pandas as pd

def extract_high_coverage_single_copy_orthogroups(orthogroups_file, output_file, coverage_threshold=0.95):
    """
    从Orthogroups.tsv中提取覆盖率达到指定阈值(默认95%)且为单拷贝的直系同源组
    
    参数:
        orthogroups_file: Orthogroups.tsv文件路径
        output_file: 输出文件路径
        coverage_threshold: 覆盖率阈值(0-1之间)
    """
    # 读取Orthogroups.tsv文件
    df = pd.read_csv(orthogroups_file, sep='\t', index_col=0)
    
    # 计算每个直系同源组的覆盖率(非空值的比例)
    coverage = df.notnull().mean(axis=1)
    
    # 检查每个直系同源组是否为单拷贝
    def is_single_copy(row):
        # 对每个物种的基因计数，逗号分隔表示多拷贝
        for genes in row.dropna():
            if ',' in str(genes):  # 如果有逗号，表示该物种有多个拷贝
                return False
        return True
    
    # 应用单拷贝检查
    single_copy_mask = df.apply(is_single_copy, axis=1)
    
    # 筛选同时满足高覆盖率和单拷贝条件的直系同源组
    high_coverage_single_copy = coverage[(coverage >= coverage_threshold) & single_copy_mask].index
    
    # 保存结果到文件
    with open(output_file, 'w') as f:
        for og in high_coverage_single_copy:
            f.write(f"{og}\n")
    
    print(f"找到 {len(high_coverage_single_copy)} 个覆盖率达到 {coverage_threshold*100}% 的单拷贝直系同源组")
    print(f"结果已保存到 {output_file}")

if __name__ == "__main__":
    import argparse
    
    parser = argparse.ArgumentParser(description='提取覆盖率达到指定阈值且为单拷贝的直系同源组')
    parser.add_argument('-i', '--input', required=True, help='Orthogroups.tsv文件路径')
    parser.add_argument('-o', '--output', required=True, help='输出文件路径')
    parser.add_argument('-t', '--threshold', type=float, default=0.95, 
                        help='覆盖率阈值(0-1之间，默认0.95)')
    
    args = parser.parse_args()
    
    extract_high_coverage_single_copy_orthogroups(args.input, args.output, args.threshold)

该python脚本是针对终端而写，使用方式为

1	python extract_high_coverage_ogs.py -i Orthogroups.tsv -o SingleCopyOrthologues.txt -t 0.95

3.进化树的美化

进化树构建完成之后美化方式有多种，MEGA，FigTree等软件，在线网站EvolView，ITOL等等，以下介绍ITOL网站，功能全面，客制化能力强
ITOL网址为：https://itol.embl.de/
该软件 2020 年 10 月由免费软件转为收费，访问模式可以使用，但不能保存注释，必须立即导出。基本的树结构调整以及字体调整可以通过 control面板实现，添加 symbol、饼图及条形图等，需要自己准备好 dataset 文件。
首先根据之前的建树结果得到树文件，RAxML运行结果中树文件名为RAxML_bipartitionsBranchLabels.out
如果想在图中添加更多的信息，就需要自己准备好配置文件，将文件拖拽到进化树美化页面。
然后可以根据需要准备以下文件：

##color_label.txt（展示文字的颜色）:
# 指定dataset类型
TREE_COLORS
# 指定文件分隔符
SEPARATOR TAB
# 指定数据，第一列为label，
# 第二列写为"label"
# 第三列为颜色
# 第四列为字体
# 第五列为字体大小
DATA
##color_range.txt（展示文字的背景色）:
##binary.txt（额外标记）:
##barchat.txt（直方图）:
##piechat.txt（饼图）:

准备好文件之后，打开ITOL，点击annotate下upload a tree，可以选择直接复制树文件内容，也可以上传文件。完成之后可以显示最初的进化树。
进化树如果是物种树一般需要定根，否则显示的是无根树，定根只需要在相应的节点点击，选择tree structure->re-root the tree here即可。

其余在右侧control panel中的各个参数可以根据自己的需要进行选择，选择好之后只需要将其他配置文件拖入网页中即可添加其他的注释信息。
其他具体的信息可以看ITOL网站的help信息。