月度归档：2013年09月

PHYLIP

发表于由
回复

1. PHYLIP简介

PHYLIP,the Phylogeny Inference Package, is a package of programs for inferring phylogenies (evolutionary trees). It can infer phylogenies by parsimony, compatibility, distance matrix methods, and likelihood. It can also compute consensus trees, compute distances between trees, draw trees, resample data sets by bootstrapping or jackknifing, edit trees, and compute distance matrices. It can handle data that are nucleotide sequences, protein sequences, gene frequencies, restriction sites, restriction fragments, distances, discrete characters, and continuous characters.

2. PHYLIP的下载和安装

下载PHYLIP，然后解压缩，得到phylip的安装文件包，其中包含doc、exe和src三个文件夹以及一个phylip.html文件。html文件结合doc中的文件在浏览器中打开是phylip详细的说明文件。当然，官网也提供了phylip 3.6的PDF版的Manual文件。

$ wget http://evolution.gs.washington.edu/phylip/download/phylip-3.695.tar.gz
$ tar zxf phylip-3.695.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/phylip-3.695/src/
$ cp Makefile.unx Makefile
$ make install

OrthoMCL介绍

发表于由
1

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的使用

发表于由
23

分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。