「生信」同源基因分析 |
您所在的位置：网站首页 › 聚类分析软件 › 「生信」同源基因分析 |
|
目录 介绍 环境配置 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,"$output"); my %family; my %choose; while (){ 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 最后 吐槽一下简书, 为什么不能生成目录???? 以后用遇到问题我再来补充 有问题欢迎交流~ |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |