真菌基因组聚类分析

复现一篇文章的工作时，发现作者使用的真菌基因组聚类分析工具是orthofinder，正好我也在学习相关分析，故做此笔记。

用到的软件：orthofinder，trimal，fasttree

这三个软件的安装均可通过conda

conda install orthofinder

conda install trimal

conda install fasttree

通过orthofinder获得各个基因组的单拷贝直系同源基因序列

以下部分参考自作者：生信石头

原文链接

---

#在一个合适的路径，创建个人工作文件；
 mkdir all && cd test 

1.准备所需物种的蛋白文件

统一后缀为(.pep/.fasta/.fa/.faa/.fas 均是Orthofinder可以识别的后缀)

注意：如果基因组具有可变剪切转录本，需要提取最长转录本进行（TBtools可）

2.运行orthofinder主程序

#暂时退至上一层目录
cd ..
#运行主程序
orthofinder -f all/  -M msa   -a 40  
#非conda安装,主程序运行使用orthofinder.py 

主要用到的相关参数介绍：

#-a 分析所用到的线程
#-f 指定文件夹（存放我们所有物种的序列） 
#-M 推断基因树的方法 可选：msa 和 dendroblast （默认 dendroblast）
dendroblast不依赖多序列比对，基于Blast评分方法聚类的方法，更节约时间。但相对多序列比对（msa）还是准确性差一点。
#-S 序列比对的方法 可选：Diamond 和 blast (默认Diamond)
diamond相对于blast比对速度更快，准确性也有保证
#-T 建树的方法 可选：fasttree, raxml, raxml-ng, iqtree (默认fasttree)
建树的精准度/耗时 raxml > iqtree > fastree; 
如果追求更高的精准度可以使用 iqtree。
# 此处应有误，最准确应该是raxml，也是最慢的 - CJ

结果解读

产生了如上图的结果文件，

• Orthogroup_Sequences 该文件夹包含了每个同源基因集合，各物种的同源基因序列。

• Orthogroups 同源组信息的目录

Orthogroups.GeneCount.tsv #每个物种在每个同源基因集合所具有的基因数目
Orthogroups.tsv #每个物种在每个同源基因集合的基因ID
Orthogroups_UnassignedGenes.tsv #每个物种在每个同源基因集合的基因ID（包括未分配同源组的基因）
Orthogroups.txt #OrthoMCL的输出格式
Orthogroups_SingleCopyOrthologues.txt #单拷贝的同源基因集合

• Single_Copy_Orthologue_Sequences 该文件包含了单拷贝的直系同源基因核酸序列。后续需要若需要构建时间分歧进化树，使用的序列。

• MultipleSequenceAlignments 多序列比对的文件。

• WorkingDirectory 运行程序的文件夹。

• Species_Tree 物种树文件夹

---

生信石头作者，程老师在结果文件中描述的是，Species_Tree文件夹中应该有如下文件

Orthogroups_for_concatenated_alignment.txt #构建进化树所用到的同源基因集合
SpeciesTree_rooted.txt #有根物种树文件
SpeciesTree_rooted_node_labels.txt
#具有Node信息的树文件；导进查看树文件的软件即可，大致了解到物种关系。

但是很可惜，我运行后，只得到了Orthogroups_for_concatenated_alignment.txt文件

要知道，我一开始选择orthofinder就是看中了其方便快捷的建树能力，可是现在没有建树结果文件了，我查看了运行日志

发现报错出现在miniconda3的bin目录下，水平有限，实在没法解决。

万幸，我检索到了另一个老师的文章：郝永超M1racle：

单拷贝直系同源基因构建物种进化树

学习了通过单拷贝直系同源基因构建物种进化树的方法

前文中，我们运行了orthofinder程序，得到了一些结果，我们要的单拷贝直系同源基因序列在 Single_Copy_Orthologue_Sequences 文件夹里

分析并建树

使用mafft软件比对，使用seqkit对齐物种并将序列转为一行，最后用paste命令把序列连起来，类似于重测序分析时把SNP连起来建树

cd ./Single_Copy_Orthologue_Sequences
ls *fa|while read file;do mafft --auto $file > $file.aln;done
ls *aln|while read file;do seqkit sort $file|seqkit seq -w 0 > $file.format;done
paste -d " " *format > all.aln.fa

得到的all.aln.fa文件中，所有序列被串联起来了，但是每条序列的名称都特别长，郝永超M1racle老师通过编写一个python脚本修改每条序列的名称，我建议手动修改，因为我还不擅长写代码，只会用。

修改后的文件命名为：format.all.aln.fa

之后，使用trimal对不符合条件的序列进行修剪，fasttree建树

~/root/trimal/source/trimal -in format.all.aln.fa -out trimal.format.all.aln.fa -automated1  #trimal修剪
fasttree <trimal.format.all.aln.fa> tree #建树

由于我用的是自己在腾讯云通过248租的三年的服务器分析的，，运行内存只有2G，储存内存只有50G，所以跑fasttree的时候出现了out of memory报错，关于out of memory报错，CSDN和简书上有很多解决方法，这里我就不赘述了。

我用别的博主介绍的释放内存的方法，释放完了，可是还是报错out of memory，这时我意识到，是这个服务器运行内存太小了。

之后，我检索fasttree的使用方法，发现了fasttree有windows版，于是，在fasttree的官网（http://www.microbesonline.org/fasttree/）下载了fasttree的windows版，

打开windows的cmd，用windows完成建树，

把刚刚下载好的fasttree和我们的输入文件trimal.format.all.aln.fa放到以下目录：C：\Users\Administrator

输入代码

fasttree <trimal.format.all.aln.fa> tree #建树

完成建树，成功得到树文件tree

随后可通过ITOL对建树结果进行美化和编辑。

完工！