目录
- 介绍
- 环境配置
- OrthoMCL使用
- 统计聚类结果
- 最后
介绍
- OrthoMCL是一款直系同源基因聚类软件, 它不仅能得到多个物种共有的直系同源基因, 还能够分别获得不同物种特有基因家族的扩张情况
-
关于直系同源和旁系同源的区别, 我就盗一张图简单解释吧
直系/旁系关系图
- 这个软件操作起来还算简单, 官方给的
userguide还算详细, 但问题是, 运行前需要确认服务器存储等各种环境, 并且还要搞定mysql数据库, 下面就从环境配置到拿到结果一步一步说吧
环境配置
因为在之前服务器上已经配过一遍了, 怕出问题, 就不在那重配了. 上个月在华为云买了个1核1G的垃圾服务器, 今天就搞它了
OrthoMCL安装
Too simple, 只需要下载解压就可以了
- 下载
$wget http://orthomcl.org/common/downloads/software/v2.0/orthomclSoftware-v2.0.9.tar.gz
- 解压
$tar zxvf orthomclSoftware-v2.0.9.tar.gz
- 添加到环境
$vi ~/.bash_profile
export PATH=/root/software/orthomclSoftware-v2.0.9/bin:$PATH
环境配置
软件运行需要UNIX系统, BLAST工具, Oracle/MySql数据库, Perl, MCL, 推荐硬件配置为4G内存+100G存储
-
BLAST安装
$conda install blast
强烈推荐使用Anaconda来下载并管理软件, 方便又条理
除了常用的BLASTN/P/X, 还有好多其他的啊, 以后用的时候再看吧...
-
MySql安装及配置
在 orthomclSoftware-v2.0.9/doc/OrthoMCLEngine/Main/mysqlInstallGuide.txt 文件中有非常详细的介绍
- 安装
$yum install -y mysql-server mysql mysql-devel
$service mysqld start #开启服务
$mysqladmin -u root password '******' #创建管理员账号和密码
$service mysqld restart
$mysql -u root -p #检查登陆
大多数做生信的服务器都会安装此数据库, 大家只需找管理员创建个人账户即 可, 如果服务器之前没有安装, 非超级用户安装起来会麻烦一些
- 利用
OrthoMCL提供的config文件进行编译
$mysql -u root -p #登陆数据库超级用户
mysql> CREATE DATABASE orthomcl; #创建database, 我的路径为 /var/lib/mysql
mysql> GRANT SELECT,INSERT,UPDATE,DELETE,CREATE VIEW,CREATE, INDEX, DROP on orthomcl.* TO orthomcl@localhost; #新建跑OrthoMCL的账号
mysql> set password for orthomcl@localhost = password('yourpassword'); #设置密码
$cd /home/scr/02_software/orthomclSoftware-v2.0.9/doc/OrthoMCLEngine/Main
$cp mysql.cnf my.cnf
$vi my.cnf
#只留下以下部分, 其他全部注释掉
[client]
[mysqld]
myisam_sort_buffer_size=4G
myisam_max_sort_file_size=200G
read_buffer_size=2G
mysql --defaults-file=my.cnf -u orthomcl -p #登陆成功即配置完成
-
Perl
没有perl的只能从安装perl开始了,
- 检查是否有
DBI and DBD::mysql modules是否安装
$perl -MDBI -e 1
$perl -MDBD::mysql -e 1
- 哈哈, 让管理员给你装吧...
perl -MCPAN -e shell
cpan> o conf makepl_arg "mysql_config=/path_to_your_mysql_dir/bin/mysql_config"
cpan> install Data::Dumper
cpan> install DBI
cpan> force install DBD::mysql
卧槽, 华为云给配的这个系统还挺好的, perl的的两个modules竟然都有, 哈哈
安装说明文档(上面提到过)中也有 non-root 用户的安装说明, 很繁琐, 反正我也没用, 不说了, 嘻嘻~
-
MCL安装
$conda install mcl
OrthoMCL使用
整个过程需要13个步骤, 下面一一介绍
-
Install and configure the relational database
通过Oracle/MySQL的userguide进行安装和配置, 上面已经做了
-
-
install mcl
推荐的是去官网下载, 我是用conda装的
-
-
install and configure OrthoMCL programs
install就不说了, 这userguide是真的啰嗦, 能读到这部分, 还能不下载&解压软件??
-
$mkdir orthomcl #创建自己的工作目录
$cd orthomcl
$cp ~/02_software/orthomclSoftware-v2.0.9/doc/OrthoMCLEngine/Main/orthomcl.config.template .out_01/00.orthomcl.config
$vi 00.orthomcl.config
# this config assumes a mysql database named 'orthomcl'. adjust according
# to your situation.
dbVendor=mysql
dbConnectString=dbi:mysql:orthomcl:localhost:3307 #设置你使用的数据库和hostname及其使用端口，默认是3307;
dbLogin=orthomcl
dbPassword=5201314
similarSequencesTable=SimilarSequences_new #以下五项可以修改的
orthologTable=Ortholog_new
inParalogTable=InParalog_new
coOrthologTable=CoOrtholog_new
interTaxonMatchView=InterTaxonMatch_new
percentMatchCutoff=50
evalueExponentCutoff=-5
oracleIndexTblSpc=NONE
orthomclInstallSchema
#将上一步设置的模型提交给 database
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema out_01/00.orthomcl.config out_01/install_tables.log
orthomclAdjustFasta
#处理 fasta 格式文件
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclAdjustFasta pdel out_02/pdel.pep.fa 1
# argv[1] 表示修改后每条序列的开头名称及文件名, argv[2]表示取原始序列名称的第一部分作为名称, 两个名称之间用'|'连接
# 将要做同源分析的物种逐个处理
orthomclFilterFasta
#序列筛选(The filter is based on length and percent stop codons)
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclFilterFasta out_02 10 20
# argv[1] 为存放上一步结果的目录, argv[2] 表示蛋白序列最小允许长度, 官方建议为10, argv[3] 表示终止密码子最大比例, 官方建议为20
All-v-all BLAST
#blastp, 最好时间, 可以拆分一下再比对
$makeblastdb -in out_03/goodProteins.fasta -dbtype prot -out out_04/orthomcl
$blastp -query out_03/goodProteins.fa -out out_04/orthomcl_blastp.out -db out_04/orthomcl -evalue 1e-5 -num_threads 5
orthomclBlastParser
#处理比对结果, 用于提交给orthomcl database
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclBlastParser out_04/orthomcl_blastp.out out_02 >> out05/similarSequences.txt
#要根据结果文件大小修改my.cnf中参数, 太累了,不想写了, 参考mysqlConfigurationGuide.txt改吧
orthomclLoadBlast
#提交给orthomcl database
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclLoadBlast out_01/00.orthomcl.config out_05/ilarSequences.txt
orthomclPairs
#这一步是主要的计算环节, 用于找到配对的蛋白
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclPairs out_01/00.orthomcl.config out_07/orthomcl_pairs.log cleanup=no
#第二次运行时往往会出错, 要把之前的database删掉, 重跑上边的命令
orthomclDumpPairsFiles
#获得ortholog, coortholog, inparalog文件
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclDumpPairsFiles out_01/00.orthomcl.config
mcl
#马尔科夫模型聚类算法软件
$mcl mclInput --abc -I 1.5 -o out_09/mclOutput
orthomclMclToGroups
#输出聚类之后的结果文件
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclMclToGroups all_ 1 < out_09/mclOutput > out10/all.txt
统计聚类结果
-
perl脚本
#!/usr/bin/perl -w
use strict;
#用于基因家族的鉴定
my $input=shift;
my $output=shift;
open (I,"<$input");
open (O,">$output");
my %family;
my %choose;
while (<I>){
chomp;
my $line=$_;
my @inf=split/\s+/,$line;
my $family=shift @inf;
my %unique;
foreach my $gene (@inf){
$gene=~/^(\w+)\|/;
my $species=$1;
$family{$species}{$family}{$gene}++;
$unique{$species}++;
}
my @species=keys %unique;
my $speciesnumber=scalar(@species);
if ($speciesnumber==1){
$choose{$species[0]}{family}++;
$choose{$species[0]}{genes}+=$unique{$species[0]};
}
}
print O "species\tgene number\tfamily number\n";
my @species=keys %family;
foreach my $species (@species){
my $totalgene=0;
my @familyid=keys %{$family{$species}};
my $familynumber=scalar(@familyid);
foreach my $familyid (@familyid){
my @genes=keys %{$family{$species}{$familyid}};
my $genenumber=scalar(@genes);
$totalgene+=$genenumber;
}
print O "$species\t$totalgene\t$familynumber\n";
}
print O "species\tuniuqe family\tunique genes\n";
my @unique=keys %choose;
foreach my $unique (@unique){
print O "$unique\t$choose{$unique}{family}\t$choose{$unique}{genes}\n";
}
close I;
close O;
- 结果展示
species gene number family number
Rcom 20608 14905
Grai 31538 15160
Egra 28651 13941
Peup 47929 17082
Atha 23426 13234
Swil 30699 17085
Ptri 35629 21303
Vvin 19252 13044
pdel 34566 20516
species uniuqe family unique genes
Grai 799 2576
Egra 842 3296
Rcom 792 2354
Atha 787 3024
Swil 223 512
Peup 375 975
Ptri 291 689
pdel 138 309
Vvin 665 1915
最后
- 吐槽一下简书, 为什么不能生成目录????
- 以后用遇到问题我再来补充
- 有问题欢迎交流~
网友评论