目錄

• 介紹
• 環(huán)境配置
• OrthoMCL使用
• 統(tǒng)計(jì)聚類結(jié)果
• 最后

介紹

• OrthoMCL是一款直系同源基因聚類軟件, 它不僅能得到多個(gè)物種共有的直系同源基因, 還能夠分別獲得不同物種特有基因家族的擴(kuò)張情況

• 關(guān)于直系同源和旁系同源的區(qū)別, 我就盜一張圖簡(jiǎn)單解釋吧

  
  
  直系/旁系關(guān)系圖
• 這個(gè)軟件操作起來(lái)還算簡(jiǎn)單, 官方給的userguide還算詳細(xì), 但問(wèn)題是, 運(yùn)行前需要確認(rèn)服務(wù)器存儲(chǔ)等各種環(huán)境, 并且還要搞定mysql數(shù)據(jù)庫(kù), 下面就從環(huán)境配置到拿到結(jié)果一步一步說(shuō)吧

環(huán)境配置

因?yàn)樵谥胺?wù)器上已經(jīng)配過(guò)一遍了, 怕出問(wèn)題, 就不在那重配了. 上個(gè)月在華為云買了個(gè)1核1G的垃圾服務(wù)器, 今天就搞它了

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

• 添加到環(huán)境

$vi ~/.bash_profile
export PATH=/root/software/orthomclSoftware-v2.0.9/bin:$PATH

環(huán)境配置

軟件運(yùn)行需要UNIX系統(tǒng), BLAST工具, Oracle/MySql數(shù)據(jù)庫(kù), Perl, MCL, 推薦硬件配置為4G內(nèi)存+100G存儲(chǔ)

• BLAST安裝

$conda install blast

強(qiáng)烈推薦使用Anaconda來(lái)下載并管理軟件, 方便又條理
除了常用的BLASTN/P/X, 還有好多其他的啊, 以后用的時(shí)候再看吧...

• MySql安裝及配置
  在 orthomclSoftware-v2.0.9/doc/OrthoMCLEngine/Main/mysqlInstallGuide.txt 文件中有非常詳細(xì)的介紹

1. 安裝

$yum install -y mysql-server mysql mysql-devel
$service mysqld start   #開啟服務(wù)
$mysqladmin -u root password '******' #創(chuàng)建管理員賬號(hào)和密碼
$service mysqld restart
$mysql -u root -p  #檢查登陸

大多數(shù)做生信的服務(wù)器都會(huì)安裝此數(shù)據(jù)庫(kù), 大家只需找管理員創(chuàng)建個(gè)人賬戶即 可, 如果服務(wù)器之前沒(méi)有安裝, 非超級(jí)用戶安裝起來(lái)會(huì)麻煩一些

2. 利用 OrthoMCL 提供的 config 文件進(jìn)行編譯

$mysql -u root -p   #登陸數(shù)據(jù)庫(kù)超級(jí)用戶

mysql> CREATE DATABASE orthomcl;   #創(chuàng)建database, 我的路徑為 /var/lib/mysql
mysql> GRANT SELECT,INSERT,UPDATE,DELETE,CREATE VIEW,CREATE, INDEX, DROP on orthomcl.* TO orthomcl@localhost;  #新建跑OrthoMCL的賬號(hào)
mysql> set password for orthomcl@localhost = password('yourpassword'); #設(shè)置密碼

$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
  沒(méi)有perl的只能從安裝perl開始了,

1. 檢查是否有 DBI and DBD::mysql modules 是否安裝

$perl -MDBI -e 1
$perl -MDBD::mysql -e 1

2. 哈哈, 讓管理員給你裝吧...

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

臥槽, 華為云給配的這個(gè)系統(tǒng)還挺好的, perl的的兩個(gè)modules竟然都有, 哈哈
安裝說(shuō)明文檔(上面提到過(guò))中也有 non-root 用戶的安裝說(shuō)明, 很繁瑣, 反正我也沒(méi)用, 不說(shuō)了, 嘻嘻~

• MCL安裝

$conda install mcl

OrthoMCL使用

整個(gè)過(guò)程需要13個(gè)步驟, 下面一一介紹

• 1. Install and configure the relational database
     通過(guò) Oracle/MySQL 的 userguide 進(jìn)行安裝和配置, 上面已經(jīng)做了
• 2. install mcl
     推薦的是去官網(wǎng)下載, 我是用conda裝的
• 3. install and configure OrthoMCL programs
install就不說(shuō)了, 這 userguide 是真的啰嗦, 能讀到這部分, 還能不下載&解壓軟件??

$mkdir orthomcl  #創(chuàng)建自己的工作目錄
$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 #設(shè)置你使用的數(shù)據(jù)庫(kù)和hostname及其使用端口，默認(rèn)是3307;
dbLogin=orthomcl
dbPassword=5201314
similarSequencesTable=SimilarSequences_new   #以下五項(xiàng)可以修改的
orthologTable=Ortholog_new 
inParalogTable=InParalog_new
coOrthologTable=CoOrtholog_new
interTaxonMatchView=InterTaxonMatch_new
percentMatchCutoff=50
evalueExponentCutoff=-5
oracleIndexTblSpc=NONE

• 4. orthomclInstallSchema

#將上一步設(shè)置的模型提交給 database 
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclInstallSchema out_01/00.orthomcl.config out_01/install_tables.log

• 5. orthomclAdjustFasta

#處理 fasta 格式文件
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclAdjustFasta pdel out_02/pdel.pep.fa 1  
# argv[1] 表示修改后每條序列的開頭名稱及文件名, argv[2]表示取原始序列名稱的第一部分作為名稱, 兩個(gè)名稱之間用'|'連接
# 將要做同源分析的物種逐個(gè)處理

• 6. 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] 為存放上一步結(jié)果的目錄, argv[2] 表示蛋白序列最小允許長(zhǎng)度, 官方建議為10, argv[3] 表示終止密碼子最大比例, 官方建議為20

• 7. All-v-all BLAST

#blastp, 最好時(shí)間, 可以拆分一下再比對(duì)
$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

• 8. orthomclBlastParser

#處理比對(duì)結(jié)果, 用于提交給orthomcl database
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclBlastParser out_04/orthomcl_blastp.out out_02 >> out05/similarSequences.txt
#要根據(jù)結(jié)果文件大小修改my.cnf中參數(shù), 太累了,不想寫了, 參考mysqlConfigurationGuide.txt改吧

• 9. orthomclLoadBlast

#提交給orthomcl database
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclLoadBlast out_01/00.orthomcl.config out_05/ilarSequences.txt

• 10. orthomclPairs

#這一步是主要的計(jì)算環(huán)節(jié), 用于找到配對(duì)的蛋白
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclPairs out_01/00.orthomcl.config out_07/orthomcl_pairs.log cleanup=no
#第二次運(yùn)行時(shí)往往會(huì)出錯(cuò), 要把之前的database刪掉, 重跑上邊的命令

• 11. orthomclDumpPairsFiles

#獲得ortholog, coortholog, inparalog文件
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclDumpPairsFiles out_01/00.orthomcl.config

• 12. mcl

#馬爾科夫模型聚類算法軟件
$mcl mclInput --abc -I 1.5 -o out_09/mclOutput

• 13. orthomclMclToGroups

#輸出聚類之后的結(jié)果文件
$/home/scr/02_software/orthomclSoftware-v2.0.9/bin/orthomclMclToGroups all_ 1 < out_09/mclOutput > out10/all.txt

統(tǒng)計(jì)聚類結(jié)果

• 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;

• 結(jié)果展示

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

最后

• 吐槽一下簡(jiǎn)書, 為什么不能生成目錄????
• 以后用遇到問(wèn)題我再來(lái)補(bǔ)充
• 有問(wèn)題歡迎交流~