参考：
利用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文件夹。
。
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 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好像选不出来保守基因，可能是单拷贝同源基因太少了。（个人推测）
个人觉得这种单纯的树没啥用，后续还要加分化时间（r8s等）还有CAFE分析。
