运行环境
ubuntu9.04
bioperl.1.6.0
paml4.4
根据 perldoc Bio::Tools::Run::Phylo::PAML::Codeml 中的例子修改
use strict;
use warnings;
use Bio::Tools::Run::Phylo::PAML::Codeml;
use Bio::AlignIO;
my $alignio = Bio::AlignIO->new(-format => 'fasta',
-file => '/home/ocrack/workspace/paml/gstd1_seqfile.fa');
my $aln = $alignio->next_aln;
my $codeml = Bio::Tools::Run::Phylo::PAML::Codeml->new();
$codeml->alignment($aln);
my ($rc,$parser) = $codeml->run();
my $result = $parser->next_result;
my $MLmatrix = $result->get_MLmatrix();
print "Ka = ", $MLmatrix->[0]->[1]->{'dN'},"\n";
print "Ks = ", $MLmatrix->[0]->[1]->{'dS'},"\n";
print "Ka/Ks = ", $MLmatrix->[0]->[1]->{'omega'},"\n";
|
出现如下错误信息:
------------- EXCEPTION: Bio::Root::NotImplemented -------------
MSG: Unknown format of PAML output did not see seqtype
STACK: Error::throw
STACK: Bio::Root::Root::throw /usr/local/share/perl/5.10.0/Bio/Root/Root.pm:357
STACK: Bio::Tools::Phylo::PAML::_parse_summary /usr/local/share/perl/5.10.0/Bio/Tools/Phylo/PAML.pm:448
STACK: Bio::Tools::Phylo::PAML::next_result /usr/local/share/perl/5.10.0/Bio/Tools/Phylo/PAML.pm:257
STACK: /home/ocrack/workspace/paml/parsCodeml2.pl:14
从
bioperl抛出的出错信息可以知道问题出在模块/usr/local/share/perl/5.10.0/Bio/Tools/Phylo
/PAML.pm中的_parse_summary函数上,当此错误发生时,调用模块Bio::Bio::Root中的throw方法处理,打印出错的原
因"MSG: Unknown format of PAML output did not see seqtype",和程序的堆栈信息.
此程序主要是对paml的结果文件进行parse,先浏览一下程序,明白程序实现过程,我是直接打开调试器(我用的是ddd和eclipse自带的调试器)运行代码,判断程序的出错原因.
在
next_result入口处设置断点,单步进入,浏览next_result的代码,明白next_result是先调用_parse_summary
函数获取paml结果文件的信息.然后通过seqtype的值判断所执行的是那个paml程序(codeml,baseml,yn00等).在调用
_parse_summary处设置断点,单步进入,浏览_parse_summary函数,明白错误应该是未parse到seqtype(取值应为
CODONML,BASE等),单步执行,分析代码,终于发现原因:
第660行
if( /^Codon\s+(usage|position)/ || /Model/)
|
使得_parse_patterns处理了原本应该_parse_summary中的正则表达式$SEQTYPE的:
CODONML (in paml version 4.4, January 2010) /tmp/w4A2DlCruV/xoHSVmFt1v
所以很自然的修改办法是
if( /^Codon\s+(usage|position)/ || /Model/ || /CODONML/ || /BASEML/ || /AAML/ || /YN00/ || /CODON2AAML/)
|
运行
perl parsCodeml2.pl
结果如下:
Ka = 0.0193
Ks = 0.1526
Ka/Ks = 0.12662
大功告成
ps.附件是修改后的PAML.pm模块和序列
|
文件: | paml.tar.gz |
大小: | 14KB |
下载: | 下载 |
|
阅读(1244) | 评论(0) | 转发(0) |