本文主要是对这篇文章利用orthofinder寻找单拷贝基因构建系统发育树 - 简书 (jianshu.com)
的复现及自己的一点心得体会。
前几天老师给了一个枯草芽孢杆菌的全基因组蛋白序列，让我去建个系统发育树，刚开始整个人都是懵的，全基因组，蛋白序列，直系同源基因，单拷贝基因。。。有点陌生，去查找了资料发现还需要自己去下载芽孢杆菌属物种蛋白序列，利用这些蛋白序列查找到直系同源单拷贝序列，对这些序列进行对比后修剪再串联构建系统发育树。主要的步骤有：
1.序列下载
在NCBI的Nucleotide数据库中以“Bacillus”，“genome”为关键词进行检索，筛选序列长度为1000000-10000000bp的序列，这样能够筛选出芽孢杆菌属物种的全基因组，再下载选定的物种的蛋白序列并将这些序列上传到服务器中。
2.orthofinder的下载与使用
2.1orthofinder下载
相比较于自己去下载orthofinder的压缩包然后自己去配置，conda安装简直不要太香，但是安装时也有点小门道，不注意可能造成orthofinder毁miniconda亡（小声bb，才不是因为自己踩过这个坑）
conda create -n orthofinder
#为orthofinder创建一个小环境，后续其依赖库都安装在这个环境中
conda activate orthofinder
#激活刚创建的环境
conda install -c bioconda orthofinder
#安装orthofinder
2.2使用
最基本的使用命令
nohup orthofinder -t 16 -f /home/monkeyflower/bioworkplace/KCYB0330 -S diamond
#nohup:不挂起运行应用，即关闭终端仍然可以继续运行。-t:指定运行线程数。-f:指定工作目录，即存放序列的目录。-S:指定序列搜索程序，有blast, mmseqs, blast_gz, diamond可用，diamond是速度最快的。
结果文件比较多，但我们在这里只需要使用Single_Copy_Orthologue_Sequences中的序列
，这个文件夹中保存的就是单拷贝直系同源序列了
3.使用muscle进行多序列比对
conda create -n muscle
conda activate muscle
conda install -c bioconda muscle
muscle -align seq.fasta -output aln.fasta
上面实现了使用miniconda设置环境，安装muscle和最基础的使用命令。
但是序列数量很多，靠一个个的去敲命令比对太浪费时间，在这里使用一点shell命令缩短时间
#!/bin/bash
# 创建一个名为aligments的文件夹
mkdir -p alignments
# 遍寻~/sequence文件夹中的所有.fa文件
for file in /home/monkeyflower/bioworkplace/KCYB0330/OrthoFinder/Results_Apr03_1/sequences/*.fa; do
# 获取没有扩展名的文件名
name=$(basename "$file" .fasta)
# 使用muscle比对序列并将比对结果输出到alignments文件夹中
muscle -align "$file" -output "alignments/${name}_aligned.fasta"
done
4.使用Gblocks进行序列修剪
依然是朴实无华的conda安装
conda create -n Gblocks
conda activate Gblocks
conda install -c bioconda gblocks
上述步骤与安装muscle时一样，先创造环境，再激活环境，再安装
批量修剪序列，代码如下
#!/bin/bash
# 进入需要进行序列修剪的目录
cd /home/monkeyflower/bioworkplace/KCYB0330/OrthoFinder/Results_Apr03_1/sequences/alignments
# 循环处理每一个序列文件
for file in *.fasta
do
# 定义输出文件名
output="/home/monkeyflower/bioworkplace/KCYB0330/OrthoFinder/Results_Apr03_1/sequences/alignments/output ${file%.fasta}_1.fasta"
# 运行 Gblocks 进行序列修剪
Gblocks "$file" -t "$output" -b4=5 -b5=h -t=p -e=.2
# 输出处理完成的序列文件名
echo "Processed $file -> /home/monkeyflower/bioworkplace/KCYB0330/OrthoFinder/Results_Apr03_1/sequences/alignments/output/"
done
5.使用seqkit串联序列
经过上一步的序列修剪后，我们需要将每个物种单个的序列进行串联，构建基因矩阵。也需要使用一点代码循环完成该任务。
conda create -n seqkit
conda activate seqkit
conda install -c bioconda seqkit
#上述代码完成了对seqkit的环境建设和安装
for i in *.2;do seqkit sort $i >$i.3
seqkit seq $i.3 -w 0 > $i.3.4
done
mkdir new
mv *.4 new/
cd new
paste -d " " *.4 > all.fa
sed -i "s\ \\g" all.fa
6.使用记事本修改各物种名称
使用记事本打开串联好的序列，使用各个物种名替代基因名。
7.使用iqtree2构建系统发育树
iqtree2 -s all.fa -m MFP -T 16 -B 1000 -o sepcific_species
