OrthoMCL及orthofinder 两种软件进行聚类分析

OrthoMCL介绍

1. OrthoMCL的用途

基于序列的相似性，OrthoMCL能将一组proteins（比如全基因组的proteins）归类到ortholog groups、in-paralogs groups和co-orthologs。

2. OrthoMCL-DB

OrthoMCL-DB包含了很多proteins，这些proteins来自一些已经完全测序的真核或原核生物的基因组。
OrthoMCL-DB将这些proteins进行了聚类，分成很多的ortholog groups。
2010.5.31，发布了OrthoMCL-DB第4版，包含 116,536个ortholog groups、1,270,853个proteins、88个真核生物基因组、16个古菌基因组、34个细菌基因组。
2011.5.31，发布了OrthoMCL-DB第5版，包含 124,740个ortholog groups、1,398,546 个proteins、150个基因组。
2013年末即将发布OrthoMCL-DB第6版。

3. OrthoMCL的两种使用方法

1. OrthoMCL-DB的官网已经将数据中的proteins进行了ortholog的聚类，其网站提供了一个工具，用于接收上传的基因组 proteins，再将这些proteins group到相应的ortholog groups中。官网提供的工具Assign your proteins to OrthoMCL Groups用于进行分析。

2. 如果要对多个基因组的proteomes进行聚类，则可以使用OrthoMCL单机版的软件来进行运算。其用法详见：OrthoMCL的使用。

4. OrthoMCL算法

1. 将多个proteomes转换成orthomcl兼容的FASTA文件。

2. 移除低质量的序列。

3. All-versus-All BLASTP with 1e-5 cutoff。即使用这些proteomes的protein sequences构建blast数据库，再将所有的这些序列和数据库进行BLASTP比对，取evalue小于1e-5的比对结果。

4. Filter by percent match length。计算比对结果的percent match length ( 所有hsp中比对上序列的长度之和 / 两条序列中短的那条序列的长度 )。取50%的cutoff值。

5. 寻找不同物种间potential ortholog pairs(两两物种的protein序列相互是best hits)；寻找同一物种内in-paralog pairs（相互之间是better hits，即对于2个序列之中的任意一条序列，和其in-paralog序列之间的evalue值 <= 这条序列和其它物种比对的evalue值）.

6. 根据上一步结果寻找co-ortholog pairs(pairs connected by orhthology and in-paralog,并且pairs之间的evalue值低于1e-5).

7. 对所有的pairs进行E-values的Normalization，以利于下一步MCL的计算。见下一部分内容，或参考OrthoMCL Algorithm Document。

8. 将所有的ortholog,in-paralog和co-ortholog pairs，以及它们的标准化后的weight值输入到MCL程序中，来进行聚类分群。MCL documentation

All-v-all BLAST
这一步你必须自己做BLAST，即：将你上一步得到的goodProteins.fasta进行多对多的blast，参数建议设置
-m 8 -F F -b 1000 -v 1000 -a 2
EXAMPLE:blastall -p blastp -i goodProteins.fasta -d goodProteins.fasta -m 8 -F F -b 1000 -v 1000 -a 2 -o all_VS_all.out.tab
**这一步事实上为MCL提供相似矩阵 **

mcl
这一步开始对上一步给出的输出文件，进行mcl操作，开始聚类
EXAMPLE: orthomclSoftware/bin/orthomclDumpPairsFile my_orthomcl_dir/orthomcl.config
输出文件为mclOutput文件
EXAMPLE： mcl my_orthomcl_dir/mclInput --abc -I 1.5 -o my_orthomcl_dir/mclOutput
这里比较重要的参数是-I 具体看mcl文档

第7步Blast，是整个过程中最关键的一步。有以下2点需要注意：

1. 数据库中的蛋白质序列数量：在OrthoMCL DB中选取和要分析的物种亲缘关系较近的几个物种的基因组，或下载其它公布的基因组，加上要分析的物种的基因组；使用这些基因组总体的蛋白质序列来构建 Blast数据库。如果只是使用要分析的物种的蛋白质序列建数据库，则inparalogs文件中成对的序列实际上是paralogs，数目比真正的 inparalogs要多很多。使用所有的OrthoMCL DB中的序列，第5版含150个基因组，信息量太大，不使用几百个核的超算或计算机集群去运行，是很不现实的。
   2. 对数据库中所有的蛋白质序列来使用blast比对到该数据库中得到结果。如果只是对要分析的物种进行Blast，则只能得到inparalogs的信息，而没有orthologs和coorthologs。

5 分析过程

5.1 输入文件格式转化

orthomcl的输入文件为fasta格式的基因或蛋白序列，fasta文件的序列名称要求以样品名开头之后接’|’分隔，之后接每个样品的序列名（如例1），而且样品名和序列名不能有重复。

命令：orthomclAdjustFasta程序，将fasta文件转换出兼容orthomcl的fasta文件使用命令： （1）orthomclAdjustFasta A(B,ref) X1(X2,X3).fa 1，结果输出为A(B，ref).fasta。（单个跑完再合并）。本文生成样品A，B和参考序列ref为例，在compliantFasta文件夹中的 序列文件名分别为：A.fa,b.fa,ref.fa。

例1：

A|gene1

ASSRKSKWQFMGARDAGAKDELRQVYGVSERTESDGAANLIHKLRAINYTLAELGQWCAYKVGQSFLSAL

B|contig1

KDELRQVYGVSERTESD

5.2 输入文件合并过滤

使用命令：orthomclFilterFasta compliantFasta/ 10 20。允许的最短的protein长度是10，stop codons最大比例为20%；生成了两个文件（2）goodProteins.fasta和poorProteins.fasta两个文

5.3 全序列比对

将上一步的goodProteins.fasta序列进行自身的多序列比对，比对使用软件为blast+，输出结果为all.m8.anno。文件太大可以拆分比对，最后合并
/share/nas2/genome/bin//blastall -b 1000 -v 1000 -a 2 -p blastp -e 1e-5 -F F -d goodProteins.fasta -i goodProteins.fasta.div1/goodProteins.fasta.f2.106.seq -o /goodProteins.fasta.div1/goodProteins.fasta.f2.106.seq.blast -m 8

cat goodProteins.fasta.f2.*.seq.blast >（3）all_VS_all.out.tab 还可以去除重复（一列，二列）最后获得 （4）all_VS_all.result

5.4 导入比对结果
将比对结果导入mysql数据库，包含以下几个步骤：
A. 将比对结果转化为规定格式，命名为similarSequences.txt，命令为：（5）orthomclBlastParser all_VS_all.result seq > similarSequences.txt
B. 将similarSequences.txt导入到数据库中，命令为：orthomclLoadBlast orthomcl.config.template similarSequences.txt

5.5 寻找paired蛋白
输入为数据库中的表SimilarSequences，和数据库的空表 InParalog, Ortholog, CoOrtholog tables；输出为对这些空表的操作，命令为：orthomclPairs orthomcl.config.template orthomcl_pairs.log cleanup=no。

5.6 将数据从mysql导出

生成（6）mcllnput文件和pairs目录。这个目录包含三个文件:
ortholog.txt, coortholog.txt, inparalog.txt。
每一个文件有三列: proteinA, protein B, their normalized score (See the Orthomcl Algorithm Document)。

命令为：orthomclDumpPairsFiles orthomcl.config.template。

5.7 使用mcl对paired蛋白聚类
命令为：mcl mclInput --abc -I 1.5 -o （7） mclOutput。

5.8 对结果编号
命令为：orthomclMclToGroups gf 1 < mclOutput > （8）groups.txt。家族名为gf_1，gf_2，gf_3...，格式如图2 。

论到直系同源基因分析的时候，大部分教程都是介绍OrthoMCL，这是2003年发表的一个工具，目前的引用次数已经达到了3000多，但这个软 件似乎在2013年之后就不在更新，而且安装时还需要用到MySQL（GitHub上有人尝试从MySQL转到sqlite）。

而OrthoFinder则是2015年出现的软件，目前已有400多引用。该软件持续更新，安装更加友好，因此我决定使用它来做直系同源基因的相关分析。

OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy提到，它的优点就是比其他的直系同源基因组的推断软件准确，并且速度还快。

此外他还能分析所提供物种的系统发育树，将基因树中的基因重复事件映射到物种树的分支上，还提供了一些比较基因组学中的统计结果。

OrthoFinder的分析过程

OrthoFinder的分析过程分为如下几步:

1. BLAST all-vs-all搜索。使用BLASTP以evalue=10e-3进行搜索，寻找潜在的同源基因。(除了BLAST, 还可以选择DIAMOND和MMSeq2)
2. 基于基因长度和系统发育距离对BLAST bit得分进行标准化。
3. 使用RBNHs确定同源组序列性相似度的阈值
4. 构建直系同源组图(orthogroup graph)，用作MCL的输入
5. 使用MCL对基因进行聚类，划分直系同源组