以下是运行codelml计算branch model的perl代码
|
use strict;
use warnings;
use Bio::Tools::Run::Phylo::PAML::Codeml;
use Bio::AlignIO;
use Data::Dumper;
use Bio::TreeIO;
# my ( $treefile, $setModel ) = @ARGV; #脚本的第一个参数是树文件名,第二个是控制omega的模型
#print "omega\tlnL\n";
my $alignio = Bio::AlignIO->new(
-format => 'fasta',
-file => 'seqfile.fa',
-dir => '.'
);
my $aln = $alignio->next_aln;
my $treeio = Bio::TreeIO->new(
'-format' => 'newick',
'-file' => 'treeH2.txt'
);
my $tree = $treeio->next_tree;
my $codeml = Bio::Tools::Run::Phylo::PAML::Codeml->new();
$codeml->alignment($aln);
$codeml->tree($tree);
#$codeml->save_tempfiles(1);
$codeml->set_parameter( "runmode", 0 );
$codeml->set_parameter( "model", 2 );
my ( $rc, $parser ) = $codeml->run();
my $result = $parser->next_result;
foreach my $model ( $result->get_NSSite_results ) {
print $model->{'likelihood'},"\n";
}
;
|
以下是树文件
|
(X02152Hom,U07178Sus,(M22585rab,((NM017025Rat,U13687Mus),(((AF070995C #1,(X04752Mus #1,U07177Rat #1)#1)#1,(U95378Sus #1,U13680Hom #1)#1)#1,(X53828OG1,U28410OG2)))))
|
出现如下错误:
Can't call method "get_NSSite_results" on an undefined value at /media/Study/evolution/paml/ex3/ex3.pl line 33, <GEN3> line 256.
at /media/Study/evolution/paml/ex3/ex3.pl line 33
单步执行调试发现,codelml.pm模块先将在tmp目录下生成,配置文件mlc,树文件,序列文件,然后运行codeml,一番调试发现原来问题出在树文件上,bioperl生成的树文件变成了这样
|
(X02152Hom,U07178Sus,(M22585rab,((NM017025Rat,U13687Mus),((("AF070995C #1",("X04752Mus #1","U07177Rat #1")#1)#1,("U95378Sus #1","U13680Hom #1")#1)#1,(X53828OG1,U28410OG2)))))
|
待bioperl生成配置文件,序列文件,树文件完毕后,暂停程序,从终端切换到这个tmp目录,运行codeml,得到错误如下:
sizeof(size_t) = 4 bytes
CODONML in paml version 4.4, January 2010
1 ambiguous codons are seen in the data:
---
528 bytes for distance
313296 bytes for conP
0 bytes for fhK
5000000 bytes for space
Seq #6 (AF070995C) is missing in the tree
原来是生成的树文件多了"(冒号)导致codeml不能识别树文件.打开vim去掉"(冒号),运行ok.
原来如此,定位到codeml.pm模块生成树文件处,488行处
|
487 $treeout->close();
488 close($temptreeFH);
|
添加如下代码
|
rename $temptreefile,$temptreefile.".bak";
open MYTREEFHOLD,$temptreefile.".bak" or die "can't open $temptreefile";
open MYTREEFHNEW,">$temptreefile" or die "can't open $temptreefile";
while (<MYTREEFHOLD>) {
s/"//g;
print MYTREEFHNEW $_;
}
close MYTREEFHOLD;
close MYTREEFHNEW;
|
再次运行脚本 ok,结果如下:
-5985.635237
附件是运行的序列,树,配置文件,和修改后的codeml.pm
|
|
文件:
|
codelml.tar.gz
|
|
大小:
|
13KB
|
|
下载:
|
下载
|
|
