参考：
利用orthofinder寻找单拷贝基因构建系统发育树 - 简书 (jianshu.com)

创建conda环境，我个人命名为orthofinder，在该环境下用conda就可安装orthofinder，muscle，mafft，seqkit，phyml软件。
其他软件prottest,RAxML可以官网下载，安装就看README或者INSTALL。比如prottest，最后还要把phyml软链到bin下，比如RAxML，make命令选一个跟自己系统匹配的运行一个就够了，不用运行全部四个。

1 Othofinder运行

首先激活环境，然后可以后台运行以下命令

orthofinder -f your/pep/fa/file/

我没加线程参数，所以十个物种的大概运行2-4个小时。生成文件中会有 Single_Copy_Orthologue_Sequences文件夹。

image.png

。

2.对生成的单拷贝同源序列进行批量处理，比对。

在Single_Copy_Orthologue_Sequences文件夹下运行

$ vi 01.muscle.sh
#!bin/bash
for i in *.fa
do muscle -in $i -out $i.1 #新版本的muscle命令变了。
#或者do muscle -align $i -output $i.1
#我至今不知道为什么第一次运行旧命令运行成功了
done

运行后多出一些后缀为.1的文件。*.1文件是比对好的序列文件。

3.提取保守序列

$ vi 02.gblock.sh
conda activate orthofinder
cd PATH/TO/Single_Copy_Orthologue_Sequences
for i in *.1
    do Gblocks $i -b4=5 -b5=h -t=p -e=.2
done

会生成一堆后缀是*.1.2的文件。

4.把序列排序，合并

$ vi 03.seqkit.sh
conda activate orthofinder
cd /jdfsbjcas1/ST_BJ/P21H28400N0232/chengxin2/02.data/10.all/pep/primary_transcripts/OrthoFinder/Results_May20/Single_Copy_Orthologue_Sequences
for i in *.2
do seqkit sort $i >$i.3
seqkit seq $i.3 -w 0 > $i.3.4
done

这里参考教程新建了文件夹new 把.4后缀的文件移到文件夹里然后

$ mkdir new
$ mv *.4 new/
$ cd new
$ paste -d " " *.4 > all.fa
$ sed -i "s\ \\g" all.fa

得到all.fa文件。

5.处理数据，将fa文件转为phy格式

vi 04.fa2phy.py
#!usr/bin/python
import re
with open('all.fa', 'r') as fin:
    sequences = [(m.group(1), ''.join(m.group(2).split()))
    for m in re.finditer(r'(?m)^>([^ \n]+)[^\n]*([^>]*)', fin.read())]
with open('all.phy', 'w') as fout:
    fout.write('%d %d\n' % (len(sequences), len(sequences[0][1])))
    for item in sequences:
        fout.write('%-20s %s\n' % item)

$ java -jar /your/software/to/prottest3-master/dist/prottest-3.4.2.jar -i all.phy -all-distributions -F -AIC -BIC -tc 0.5 -threads 24 -o prottest.out

得到蛋白模型。

6.RAxML构树

$ vi RAxML.sh
/your/path/to/01.Software/standard-RAxML-master/raxmlHPC-PTHREADS-SSE3 -T 16 -f a -x 123 -p 123 -N 1000 -m PROTGAMMAIJTTF -k -O -o Ensete_ventricosum,Musa_a,Musa_b \
       -n all.rename.tre -s all_rename.fa

qsub提交任务就可。十个物种大概运行1小时。
这里-o选项是定根，我的外群是芭蕉。如果不定外群，结果会很奇怪。
还有就是。-s参数后接的fa文件，其实是把每个物种的单拷贝同源序列合并成超长一条。最好把ID改成相应的物种。不然ID太长也会报错。
然后我用的besttree复制到figtree软件或者iTOL，就可以出图了。

注：我只选了同一“目”下的一些物种，我要研究的五种物种的基因组，加上近缘物种选了五种，一共是十种。因为很近，所以gblocks可以运行。因为我想在更多的物种比如单子叶和双子叶选了13种的时候gblocks好像选不出来保守基因，可能是单拷贝同源基因太少了。（个人推测）