参考了众多文章并结合实际操作后的感想。
参考: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 comliantFasta/ 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 #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工作目录
手动把blastresult文件中的gnl|替换为空,gredit中的ctrl+h命令
$ 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/\t(\w+)(\|.*)orthomcl/\t$1$2$1/‘ similarSequences.txt #将subject_taxon修改为正确的分类名
$ perl -p -i -e ‘s/0\t0/1\t-181/‘ similarSequences.txt #修改evalue_mant, evalue_exp,将evalue为0修改为1e-181,为在后续步骤寻找pairwise relationships
9.将similarSequences.txt文件中的相似性序列载入mysql数据库中[当前路径为my_orthomcl_dir]
$ orthomclLoadBlast orthomcl.config similiarSequences.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即为最终结果文件,可对其进行各种数据操作