Orthomcl的详细使用

参考了众多文章并结合实际操作后的感想。

参考：http://www.plob.org/2013/09/18/6174.html

参考：http://www.plob.org/2012/06/12/2207.html

Orthomcl需要的linux详细配置，简单叙述。

　　系统：unix

　　BLAST

　　Database：oracle 或 mysql

　　内存4G，硬盘100G

　　perl: DBI module 和 DBD-mysql

　　MCL程序

详细操作过程及需要注意的点：

1.数据库安装和配置

　　我用的mysql，为避免mysql运行出问题，最好修改mysql配置文件：

　　1.设置myisam_sort_buffer_size为可用内存的一半

　　2.设置myisam_max_sort_file_size为orthomclBlastParser程序生成文件similarSequences.txt的5倍大小

　　3.设置read_buffer_size的值为???，或者尽量大

2.安装MCL

　　下载地址在 http://www.micans.org/mcl/src/mcl-latest.tar.gz，自动获取最新版

　　$ ./configure

　　$ make

　　$ make install

　　#注意出现make[] Nothing to be done for '***'

　　　　　　　make[] leaving directory '/*/*/***'

　　不是make错误！make时最好用root权限

3.安装配置Orthomcl software

　　下载地址 http://orthomcl.org/common/downloads/software/v2.0/ 我下的是orthomclSoftware-v2.0.9.tar.gz

　　下载完成后解压文件 $ tar -xf orthomclSoftware-v2.0.9.tar.gz，解压完成后包括bin config doc lib 四个文件夹

　　把bin目录加到环境变量里，方便以后操作

　　$ vim .bashrc #在~目录下执行此操作

　　在最后一行加 export PATH=$PATH:/home/abc/.../orthomclSoftware/bin

　　$source .bashrc

　　在Orthomcl software主文件夹下创建文件夹my_orthomcl_dir作为工作目录，把/doc/OrthoMCLEngine/Main/orthomcl.config.template文件复制到my_orthomcl_dir下，命令如下：[当前路径在my_orthomcl_dir]

　　$ mkdir my_orthomcl_dir

　　$ cd my_orthomcl_dir

　　$ cp ../doc/OrthoMCLEngine/Main/orthomcl.config.template orthomcl.config #重命名为orthomcl.config

　　打开orthomcl.config：

　　##

　　dbVendor=mysql　　#使用的数据库为mysql，若为oracle，则改成oracle

　　dbConnectString=dbi:mysql:orthomcl　　#连接到mysql里的orthomcl数据库

　　dbLogin=test　　#数据库的用户名

　　dbPassword=123　　#与用户名相对应的密码

　　orthologTable=Ortholog

　　inParalogTable=InParalog

　　coOrthologTable=CoOrtholog

　　interTaxonMatchView=InterTaxonMatch

　　percentMatchCutoff=50

　　evalueExponentCutoff=-5

　　oracleIndexTblSpc=NONE

4.在mysql中创建数据库，并创建表

　　$ mysql -u root -p　　#先用root登录创建名为orthomcl的数据库（如果test用户有权限创建数据库的话就用test用户登录创建）

　　mysql> create database orthomcl;

　　mysql> grant all on orthomcl.* to 'test'@'%';　　#赋予test用户操作orthomcl数据库的所有权限，test处为用户名，%代表任何主机

　　$ orthomclInstallSchema orthomcl.config mysql.log species　　#执行orthomclInstallSchema命令，按照orthomcl.config的配置，在数据库中创建表，mysql.log记录日志文件（选） species为每个表名后加物种名（选）

5.格式化orthomcl输入文件

　　orthomcl的输入文件要求为fasta格式，且需满足一定的要求，格式要求为如下：

　　　　>taxoncode|unique_protein_id　　#taxoncode为物种代码一般为3-4个字母；unique_protein_id为蛋白id，taxoncode和蛋白id之间用|隔开

　　 　　 MFAXGETHFD..........

　　　如：

　　　　>Dha|CAG25565

　　　　  MKFGHA...........

　　使用orthomclAdjustFasta程序可以把从NCBI上下载下来的序列文件转换成orthomcl的标准格式，转换格式前先在my_orthomcl_dir目录下创建名为compliantFasta的文件夹，命令如下：[当前路径为my_orthomcl_dir]

　　$ mkdir compliantFasta

　　$ cd compliantFasta

　　$ orthomclAdjustFasta hsa ../Homo_sapiens.NCBI36.53.pep.all.fa 1　　#hsa为物种代码;***.fa为序列文件，存放在my_orthomcl_dir目录下;1代表在id前加物种名和|

　　执行完上述命令后，产生的文件为hsa.fasta存放在compliantFasta目录下，格式如下：

　　　　>hsa|TIG25566

　　　　  MHSJKXG..........

　　compliantFasta文件夹下存放各个物种的蛋白组，如Hsa.fasta Dha.fasta Ali.fasta Kla.fasta......

　　当然，，，自己编写的脚本也很容易完成格式转换，，，

6.过滤序列

　　使用orthomclFilterFasta命令对compliantFasta文件夹下的序列进行过滤，orthomcl的推荐规则是允许protein序列最短长度为10，stop coden占的最大比例为20%，命令会在my_orthomcl_dir目录下产生goodProteins.fasta和poorProteins.fasta，goodProteins.fasta文件中包含所有comliantFasta文件夹下经过筛选的物种蛋白组。[当前路径为my_orthomcl_dir]

　　$ orthomclFilterFasta compliantFasta/ 10 20

7.创建blast数据库，对goodProteins.fasta做blastp

　　下载最新版本blast+,安装成功根据需求下载蛋白数据库，以我自己为例，我只需要compliantFasta文件夹下的几个物种蛋白组，即只需要goodProteins.fasta作为数据库

　　把blast添加到环境变量后，执行如下命令：[当前路径为blast主目录]

　　$ mkdir db　　#创建文件夹db，作为blast数据库文件夹

　　$ cd db

　　$ makeblastdb -in goodProteins.fasta -dbtype prot -title orthomcl -parse_seqids -out orthomcl -logfile orthomcl.log　　#以goodProteins.fasta为序列文件，创建名为orthomcl的blast数据库[当前路径为db]

　　$ blastp -db orthomcl -query goodProteins.fasta -seg yes -out orthomcl.blastout -evalue 1e-5 -outfmt 7 -num_threads 24　　#goodProteins.fasta对orthomcl库做blast，产生的结果文件为orthomcl.blastout（自身对自身做blast，找出目标物种间的同源基因）

8.处理blast产生的结果（参考）[当前路径为db]

　　$ grep -P "^[^#]" orthomcl.blastout > blastresult　　#将orthomcl.blastout中的注释行去掉，产生新文件blastresult。

　　$ cp blastresult /home/abc/.../orthomclSoftware/my_orthomcl_dir　　#把blastresult复制到orthomcl工作目录下

　　$ cd /home/abc/.../orthomclSoftware/my_orthomcl_dir　　#转到orthomcl工作目录

　　$ sed -i s/"gnl|"/""/g blastresult

　　$ orthomclBlastParser blastresult compliantFasta > similarSequences.txt　　#使用orthomclBlastParser命令引入compliantFasta文件夹下文件，生成similarSequences.txt文件，找出相似性序列，从第1列到第8列分别是：query_id, subject_id, query_taxon, subject_taxon, evalue_mant, evalue_exp, percent_ident, percent_match。

　　$ perl -p -i -e 's/ (w+)(|.*)orthomcl/ $1$2$1/' similarSequences.txt　　#将subject_taxon修改为正确的分类名

　　$ perl -p -i -e 's/0 0/1 -181/' similarSequences.txt　　#修改evalue_mant, evalue_exp,将evalue为0修改为1e-181，为在后续步骤寻找pairwise relationships

9.将similarSequences.txt文件中的相似性序列载入mysql数据库中[当前路径为my_orthomcl_dir]

　　$ orthomclLoadBlast orthomcl.config similarSequences.txt　　#将数据导入数据库中

10.寻找成对蛋白质[当前路径为my_orthomcl_dir]

　　$ orthomclPairs orthomcl.config orthomcl_pairs.log cleanup=no　　#此命令对数据库中的空表进行操作

11.将数据从mysql数据库中导出[当前路径为my_orthomcl_dir]

　　$ orthomclDumpPairsFiles orthomcl.config　　#此命令会在my_orthomcl_dir下生成一个mclInput文件和一个pairs文件夹，pairs文件夹下包含coorthologs.txt和inparalogs.txt和orthologs.txt三个文件

12.使用mcl对pairs进行聚类[当前路径为my_orthomcl_dir]

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

13.提取mcl的结果，生成group.txt文件[当前路径为my_orthomcl_dir]

　　$ orthomclMclToGroups Fungi 1 < mclOutput > groups.txt　　#生成groups.txt文件，每个同源组的编号从Fungi1开始，依次递增

至此orthomcl程序运行完毕，产生的groups.txt即为最终结果文件，可对其进行各种数据操作