[原创]bioperl中解析paml结果模块的一个bug

运行环境
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
下载:	下载