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

5. pairs的evalue计算和标准化

pairs的evalue计算:pairs的两条序列相互blast后有两个evalue值,这两个值常常不相等。但是为了计算需要,于是pairs的之间的两个evalue值要进行一个计算,得到pairs weight,weight= ( -log10(evalue1) + -log10(evalue2) ) / 2 。
pairs的evalue的标准化:1. 对于in-paralog pairs,在某一个基因组中,取两条序列中任意一条序列有ortholog的in-paralog pairs为有效in-paralog pairs。若在这个基因组没有这样的pairs,则该基因组所有的in-paralog pairs都为有效的in-paalog pairs。最后得到所有基因组所有有效的in-paralog pairs。然后取这些有效in-paralog pairs的weight的平均值。最后,每个in-paralog pair的evalue标准化后的值为其weight除以average weight。 2. 对于ortholog或co-ortholog pairs则简单很多,求所有weight的平均值,然后使用各个pair的weight除以average weight,则将其标准化了。

6. 网络版的OrthoMCL的使用

OrthoMCL-DB已经对150个proteomes进行了OrthoMCL的分析,对orthologs进行了聚类。这个过程由于数据量大,因此,在好几百的CPU资源下也需要好几个星期才能做完。
在OrthoMCL-DB上上载 a set of proteins,服务器则会将所有上传的proteins比对到OrthoMCL-DB中所有的proteins上;选取evalue 1e-5和50% match的cutoff;然后将protein归类到其top hit所对应的protein的类上;如果top hit所对应的proteins没有group,则该protein归类到NO_GROUP。
然后,再对上一步cutoff掉的proteins来使用OrthoMCL-DB的in-paralog算法来创建in-paralogs pairs,然后再进行MCL的聚类。
使用该方法,最后将a set of proteins进行了同源基因的聚类,但是缺点如下:
1. 这种方法是单向最佳,根据protein比对的最佳结果去归类到已有的group中去。但是反过来,最佳比对结果对应的protein不一定和query protein是最佳的。这和OrthoMCL的算法是有出入的,所以该方法省了时间,但是结果和真正的结果是有一定差别的。
2. 只使用cutoff后剩下的proteins进行in-paralog分析,而没有进行所有query proteins之间的in-paralog分析。
3. 没有ortholog pairs和co-ortholog pairs的信息,没法进行单拷贝同源基因的提取与分析。

7. 本地OrthoMCL的使用

对指定的a set of proteomes进行同源基因分析,则使用本地的OrthoMCL进行分析。而官网不提供这种服务,因为消耗的计算机资源过大。

8. 注意事项

1. 序列都是使用protein序列,而不是nucleotide序列,是因为protein序列更精确。
2. proteomes中的序列要去除可变剪切,只留取alternative proteins中长度最长的。否则在有alternative proteins存在的情况下,则会造成pseudo-in-paralogs(即alternative proteins称为in-paralogs),给后续的分析造成麻烦。
3. paralog分为in-paralog和out-paralog。in-paralog是指同一个物种的paralogs的分化发生在物种分化之后,这样的话,代表in-paralogs之间的序列相似性比其orthologs的相似性要高;通过OrthoMCL的原理可以看出,很好得到分析。而out-paralog是指paralogs的分化发生在物种分化以前,这代表out-paralogs的序列之间的相似性比某一个物种的orthologs的相似性要低;这样是很不好分析的,因为不好定阈值,或者得到的结果不易得不到大众的认可;OrthoMCL也没进行out-paralog的分析;当然,也可以只将一个query proteome输入到orthomcl来进行分析,得到的是所有paralog分析结果,包含了out-paralog。

OrthoMCL的使用

分13步进行,如下:

1. 安装和配置数据库

Orthomcl可以使用Oracle和Mysql数据库,而在这里只介绍使用Mysql数据库。
修改配置文件/etc/my.cnf,对Mysql进行如下配置:

1. 设置myisam_sort_buffer_size的值为可用内存的一半;
2. 设置myisam_max_sort_file_size为orthomclBlastParser程序生成文件similarSequences.txt的5倍大小;
3. 软件的说明文档中设置read_buffer_size的值为???,但是设置为3个问号或1个问号,则mysql启动不了。我将其设置为2000M。

2. 安装mcl软件

mcl,即Markov Clustering algorithm,其最新的软件下载地址:http://www.micans.org/mcl/src/mcl-latest.tar.gz。下载后使用’./configure && make && make install’安装即可。

3. 安装并配置OrthoMCL软件

下载OrthoMCL软件(http://orthomcl.org/common/downloads/software/)后,解压缩后,其中包含文件夹:bin、config、doc、lib四个文件夹。
在一个工作目录中运行OrthoMCL,该目录包含数据文件和结果文件。将doc/OrthoMCLEngine/Main/orthomcl.config.template复制到工作目录中。该文件为OrthoMCL的配置文件,以使用mysql数据库为例,其中的内容如下:

dbVendor=mysql   #使用的数据库为mysql
dbConnectString=dbi:mysql:orthomcl   #使用mysql数据库中名为orthomcl的数据库
dbLogin=test    #数据库的用户名
dbPassword=123  #相应的密码
similarSequencesTable=SimilarSequences #
orthologTable=Ortholog
inParalogTable=InParalog
coOrthologTable=CoOrtholog
interTaxonMatchView=InterTaxonMatch
percentMatchCutoff=50
evalueExponentCutoff=-5
oracleIndexTblSpc=NONE

4. 安装orthomcl数据库的表

首先,进入mysql数据库,新建一个名为orthomcl的数据库;然后,使用orthomclInstallSchema命令在数据库中创建一些表,这些表的名字则是orthomcl.config.template配置文件中指定的5个名称。

$ mysql -u test -p
mysql> create database orthomcl;
$ cp /opt/biosoft/orthomclSoftware-v2.0.9/doc/OrthoMCLEngine/Main/orthomcl.config.template .
$ orthomclInstallSchema orthomcl.config.template [log species]

orthomclInstallSchema命令后面不接参数则给出帮助文档。以上命令使用orthomcl.config.template配置文件中的设置生成了数据库中相应的表,如果加入方括号中的内容,则会生成日志文件log,生成的表后缀都是species。

5. 创建orthomcl的输入文件

orthomcl的输入文件为fasta格式文件,但是fasta文件的序列名称要满足这样的要求:

>taxoncode|unique_protein_id
MHDR...
>hsa|sequence_1
MHDR...
>led|scaffold_1.1
MHDR...

序列名第一列是物种的代码,一般是3到4个字母;中间使用’|'符号隔开;第二列是蛋白质序列独一无二的id。
一般输入文件是fasta格式,其序列名由空格或’|'隔开,使用orthomclAdjustFasta程序,将fasta文件转换出兼容orthomcl的fasta文件

$ redun_remove protein.fasta > non_dun_protein.fasta
$ mkdir compliantFasta; cd compliantFasta
$ orthomclAdjustFasta led ../non_dun_protein.fasta 1

上述命令去除可变剪切的蛋白质序列;创建了文件夹compliantFasta;然后使用orthomclAdjustFasta命令选取了protein.fasta序列名的第一列作为输出的fasta文件的序列id;输出的文件为led.fasta.

6. 过滤序列

对compliantFasta文件夹中的序列进行过滤,允许的最短的protein长度是10,stop codons最大比例为20%;生成了两个文件goodProteins.fasta和poorProteins.fasta两个文件。

$ orthomclFilterFasta compliantFasta/ 10 20

compliantFasta只能过滤低质量的序列。而实际上最好还需要过滤掉可变剪切,只留取可变剪切中最长的蛋白质序列,这个需要自行解决。

7. 对goodProteins.fasta中的序列进行BLAST

下载最新版本的Blast+,和最新版本的OrthoMCL DB的protein序列,将OrthoMCL DB的protein序列加上gooProtein.fasta中的序列合到一起做成一个blast+的数据库。然后对基因组的蛋白质序列进行比对。

$ /opt/biosoft/ncbi-blast-2.2.28+/bin/makeblastdb -in orthomcl.fasta -dbtype prot -title orthomcl -parse_seqids -out orthomcl -logfile orthomcl.log
$ /opt/biosoft/ncbi-blast-2.2.28+/bin/blastp -db orthomcl -query goodProteins.fasta -seg yes -out orthomcl.blastout -evalue 1e-5 -outfmt 7 -num_threads 24

生成orthomcl的blast DB需要97秒左右;使用-outfmt 7生成带注释的表格结果,这一步需要很长时间了,取决于电脑的运算性能。我使用24个线程,每分钟运行约27.75条序列,大约7.2个小时,运行1.2万条protein序列的比对。
blast中使用了-seg yes表示使用seg程序来进行过滤,将那些影响比对结果的低复杂度区域过滤掉。blast生成的文件结果,从第1列到第12列分别是:query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, q. end, evalue, bit score。

8. 处理Blast的结果

对上一步blast的结果进行处理,从而得到序列的相似性结果,以用于导入到orthomcl数据库中。compliantFasta文件夹中包含下载下来的OrthoMCL DB的所有蛋白质数据的文件orthomcl.fasta.

$ grep -P "^[^#]" orthomcl.blastout > blastresult
$ orthomclBlastParser blastresult compliantFasta > similarSequences.txt
$ perl -p -i -e 's/\t(\w+)(\|.*)orthomcl/\t$1$2$1/' similarSequences.txt
$ perl -p -i -e 's/0\t0/1\t-181/' similarSequences.txt

第一条命令将orthomcl.blastout中的注释行去掉,生成新的文件blastresult,不然再下一个命令中会报错的。
第二条命令生成文件similarSequences.txt,从第1列到第8列分别是:query_id, subject_id, query_taxon, subject_taxon, evalue_mant, evalue_exp, percent_ident, percent_match。值得注意的是subject_taxon是orthomclBlastParser读取的在compliantFasta文件夹中fasta文件的前缀,在此结果中,这一列则全是orthomcl。
第三条命令将subject_taxon修改为正确的分类名。
第四条命令修改evalue_mant, evalue_exp,将evalue为0修改为1e-181,这在后续步骤寻找pairwise relationships时候有要求。

9. 将similarSequences.txt载入到数据库中

生成的similarSequences.txt文件大小为83M,则修改/etc/my.cnf文件,在myisam_sort_buffer_size这一行上加一行‘myisam_max_sort_file_size = 424960’。数值是83M的5倍。然后运行:

$ orthomclLoadBlast orthomcl.config.template similarSequences.txt

10. 寻找成对的蛋白质

$ orthomclPairs orthomcl.config.template orthomcl_pairs.log cleanup=no

输入为数据库中的表SimilarSequences,和数据库的空表InParalog, Ortholog, CoOrtholog tables;输出为对这些空表的操作。故配置文件中的用户要有 update/insert/truncate权限。

11. 将数据从数据库中导出

$ orthomclDumpPairsFiles orthomcl.config.template

生成了一个文件mclInput和一个文件夹pairs;文件夹中包含3个文件coorthologs.txt,inparalogs.txt,orthologs.txt。

12. 使用mcl进行对pairs进行聚类

$ mcl mclInput --abc -I 1.5 -o mclOutput

13. 对mcl的聚类结果进行编号

$ orthomclMclToGroups led 1 < mclOutput > groups.txt

对聚类结果进行编号,依次为led1,led2, led3…

注意事项

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