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

时间:2024-02-23 10:28:28

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。

 

4 分析过程

4.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

4.2 输入文件合并过滤

使用命令:orthomclFilterFasta  compliantFasta/  10  20。允许的最短的protein长度是10,stop codons最大比例为20%;生成了两个文件(2)goodProteins.fasta和poorProteins.fasta两个文

4.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

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

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

4.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。

4.7 使用mcl对paired蛋白聚类
命令为:mcl mclInput --abc -I 1.5 -o (7) mclOutput。

4.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对基因进行聚类,划分直系同源组



 

 

来源:

http://www.chenlianfu.com/?p=1839

https://orthomcl.org/orthomcl/?rm=proteomeUploadForm 

https://www.jianshu.com/p/16e0bbb2ba19

https://blog.csdn.net/ych99ych99/article/details/17436367?utm_source=blogxgwz0

https://www.plob.org/article/6174.html

https://www.cnblogs.com/renping/p/10126229.html