Chinaunix首页 | 论坛 | 博客

ocrackocrack.blog.chinaunix.net

首页　| 　博文目录　| 　关于我

ocrack

  • 博客访问： 29244
  • 博文数量： 15
  • 博客积分： 650
  • 博客等级： 上士
  • 技术积分： 160
  • 用 户 组： 普通用户
  • 注册时间： 2009-05-31 22:01
文章分类

全部博文（15）

文章存档

2011年（2）

2010年（6）

2009年（7）

我的朋友
最近访客
推荐博文
[原创]bioperl中PAML.pm 模块的一个烦人bug

分类：

2010-10-22 19:30:36

下面的程序是用于解析codeml程序计算正选择位点时产生的结果文件.

use Bio::Tools::Phylo::PAML;
use Data::Dumper;
my $paml = Bio::Tools::Phylo::PAML->new( -file => 'results.txt' );
my $result = $paml->next_result;
#print Dumper($result);

foreach my $model ( $result->get_NSSite_results ) {
    print $model->likelihood, "\n";

    # otherwise query NEB and BEB slots

print "NEB\n";
    for my $sites ( $model->get_NEB_pos_selected_sites ) {
        print join( "\t", @$sites ), "\n";
    }
print "BEB\n";
    for my $sites ( $model->get_BEB_pos_selected_sites ) {
        print join( "\t", @$sites ), "\n";
    }

}


运行结果
--------------------- WARNING ---------------------
MSG: no node could be found for 53 (no lca?)
---------------------------------------------------
Can't call method "is_Leaf" on an undefined value at /usr/local/share/perl/5.10.0/Bio/Tools/Phylo/PAML.pm line 1061, line 796.

调试很久才发现问题是PAML.pm模块的这一段

        while( my ($k,$v) = each %branch_dnds) {
         # we can probably do better by caching at some point

         my @nodes;
         for my $id ( split(/\.\./,$k ) ) {
            my @nodes_L = map { $tree->find_node(-id => $_) } @{$idlookup->{$id}};
            my $n = @nodes_L < 2 ? shift(@nodes_L) : $tree->get_lca(@nodes_L);
            if( ! $n ) {
             $self->warn("no node could be found for $id (no lca?)");
            }
            unless( $n->is_Leaf && $n->id) {
                $n->id($id);
                }
            push @nodes, $n;
         }
         my ($parent,$child) = @nodes;
         while ( my ($kk,$vv) = each %$v ) {
            $child->add_tag_value($kk,$vv);
         }
        }

先看看结果文件中的这一段
TREE #  1:  (1, 2, ((3, 5), (6, (4, (7, (8, (((20, (9, 10)), (((11, 12), (15, (13, 14))), (16, ((42, (23, (21, 22))), ((40, ((44, ((26, (24, 25)), (45, (28, 29)))), ((27, (33, (32, (30, 31)))), ((34, 35), ((36, 37), (38, 39)))))), (41, 43)))))), (19, (17, 18)))))))));   MP score: -1
lnL(ntime:  0  np:  5):  -9241.335890      +0.000000
 5.089639 0.511336 0.380850 0.052798 2.997761

tree length =   8.78314

(1: 0.014586, 2: 0.041868, ((3: 0.070966, 5: 0.011939): 0.059490, (6: 0.064122, (4: 0.092894, (7: 0.211820, (8: 0.288564, (((20: 0.225787, (9: 0.033320, 10: 0.077065): 0.195466): 0.032066, (((11: 0.022673, 12: 0.020032): 0.198520, (15: 0.139394, (13: 0.008490, 14: 0.012848): 0.244866): 0.051892): 0.014580, (16: 0.383260, ((42: 0.471532, (23: 0.238336, (21: 0.000004, 22: 0.012619): 0.187342): 0.068532): 0.045566, ((40: 0.219803, ((44: 0.184943, ((26: 0.136296, (24: 0.211558, 25: 0.106960): 0.049198): 0.091821, (45: 0.268222, (28: 0.120440, 29: 0.178342): 0.046454): 0.011009): 0.000000): 0.018553, ((27: 0.220961, (33: 0.218276, (32: 0.127441, (30: 0.077969, 31: 0.066965): 0.032972): 0.053337): 0.065516): 0.034062, ((34: 0.000004, 35: 0.016312): 0.168746, ((36: 0.027318, 37: 0.014620): 0.105226, (38: 0.085591, 39: 0.108868): 0.026650): 0.060744): 0.076203): 0.048228): 0.047638): 0.051495, (41: 0.244808, 43: 0.311371): 0.061743): 0.046768): 0.029504): 0.036780): 0.024922): 0.027657, (19: 0.221857, (17: 0.229061, 18: 0.097006): 0.045976): 0.130815): 0.036083): 0.058565): 0.216554): 0.027679): 0.012174): 0.004641);

(AJ344371: 0.014586, AJ344373: 0.041868, ((AJ344385: 0.070966, AJ344386: 0.011939): 0.059490, (AJ344415: 0.064122, (AJ344379: 0.092894, (AJ344369: 0.211820, (AJ344409: 0.288564, (((AJ344405: 0.225787, (AJ344370: 0.033320, AJ344377: 0.077065): 0.195466): 0.032066, (((AJ344382: 0.022673, AJ344383: 0.020032): 0.198520, (AJ344378: 0.139394, (AJ344376: 0.008490, AJ344398: 0.012848): 0.244866): 0.051892): 0.014580, (AJ344414: 0.383260, ((AJ344372: 0.471532, (AJ344381: 0.238336, (AJ344395: 0.000004, AJ344396: 0.012619): 0.187342): 0.068532): 0.045566, ((AJ344400: 0.219803, ((AJ344374: 0.184943, ((AJ344410: 0.136296, (AJ344384: 0.211558, AJ344391: 0.106960): 0.049198): 0.091821, (AJ344375: 0.268222, (AJ344388: 0.120440, AJ344406: 0.178342): 0.046454): 0.011009): 0.000000): 0.018553, ((AJ344408: 0.220961, (AJ344407: 0.218276, (AJ344394: 0.127441, (AJ344393: 0.077969, AJ344404: 0.066965): 0.032972): 0.053337): 0.065516): 0.034062, ((AJ344390: 0.000004, AJ344411: 0.016312): 0.168746, ((AJ344392: 0.027318, AJ344397: 0.014620): 0.105226, (AJ344399: 0.085591, AJ344401: 0.108868): 0.026650): 0.060744): 0.076203): 0.048228): 0.047638): 0.051495, (AJ344387: 0.244808, AJ344380: 0.311371): 0.061743): 0.046768): 0.029504): 0.036780): 0.024922): 0.027657, (AJ344402: 0.221857, (AJ344389: 0.229061, AJ344403: 0.097006): 0.045976): 0.130815): 0.036083): 0.058565): 0.216554): 0.027679): 0.012174): 0.004641);

Detailed output identifying parameters

kappa (ts/tv) =  5.08964


dN/dS (w) for site classes (K=3)

p:   0.51134  0.38085  0.10781
w:   0.05280  1.00000  2.99776

dN & dS for each branch

 branch          t       N       S   dN/dS      dN      dS  N*dN  S*dS

  46..1       0.015    527.4    228.6   0.7310   0.0044   0.0060    2.3    1.4
  46..2       0.042    527.4    228.6   0.7310   0.0126   0.0172    6.6    3.9
  46..47      0.005    527.4    228.6   0.7310   0.0014   0.0019    0.7    0.4
  47..48      0.059    527.4    228.6   0.7310   0.0178   0.0244    9.4    5.6
  48..3       0.071    527.4    228.6   0.7310   0.0213   0.0291   11.2    6.7
  48..5       0.012    527.4    228.6   0.7310   0.0036   0.0049    1.9    1.1
  47..49      0.012    527.4    228.6   0.7310   0.0037   0.0050    1.9    1.1
  49..6       0.064    527.4    228.6   0.7310   0.0192   0.0263   10.1    6.0
  49..50      0.028    527.4    228.6   0.7310   0.0083   0.0114    4.4    2.6
  50..4       0.093    527.4    228.6   0.7310   0.0279   0.0381   14.7    8.7
  50..51      0.217    527.4    228.6   0.7310   0.0650   0.0889   34.3   20.3
  51..7       0.212    527.4    228.6   0.7310   0.0635   0.0869   33.5   19.9
  51..52      0.059    527.4    228.6   0.7310   0.0176   0.0240    9.3    5.5
  52..8       0.289    527.4    228.6   0.7310   0.0866   0.1184   45.7   27.1
  52..53      0.036    527.4    228.6   0.7310   0.0108   0.0148    5.7    3.4
  53..54      0.028    527.4    228.6   0.7310   0.0083   0.0113    4.4    2.6
  54..55      0.032    527.4    228.6   0.7310   0.0096   0.0132    5.1    3.0
  55..20      0.226    527.4    228.6   0.7310   0.0677   0.0926   35.7   21.2
  55..56      0.195    527.4    228.6   0.7310   0.0586   0.0802   30.9   18.3
  56..9       0.033    527.4    228.6   0.7310   0.0100   0.0137    5.3    3.1
  56..10      0.077    527.4    228.6   0.7310   0.0231   0.0316   12.2    7.2
  54..57      0.025    527.4    228.6   0.7310   0.0075   0.0102    3.9    2.3
  57..58      0.015    527.4    228.6   0.7310   0.0044   0.0060    2.3    1.4
  58..59      0.199    527.4    228.6   0.7310   0.0595   0.0815   31.4   18.6
  59..11      0.023    527.4    228.6   0.7310   0.0068   0.0093    3.6    2.1
  59..12      0.020    527.4    228.6   0.7310   0.0060   0.0082    3.2    1.9
  58..60      0.052    527.4    228.6   0.7310   0.0156   0.0213    8.2    4.9
  60..15      0.139    527.4    228.6   0.7310   0.0418   0.0572   22.1   13.1
  60..61      0.245    527.4    228.6   0.7310   0.0735   0.1005   38.7   23.0
  61..13      0.008    527.4    228.6   0.7310   0.0025   0.0035    1.3    0.8
  61..14      0.013    527.4    228.6   0.7310   0.0039   0.0053    2.0    1.2
  57..62      0.037    527.4    228.6   0.7310   0.0110   0.0151    5.8    3.4
  62..16      0.383    527.4    228.6   0.7310   0.1150   0.1573   60.6   35.9
  62..63      0.030    527.4    228.6   0.7310   0.0089   0.0121    4.7    2.8
  63..64      0.046    527.4    228.6   0.7310   0.0137   0.0187    7.2    4.3
  64..42      0.472    527.4    228.6   0.7310   0.1414   0.1935   74.6   44.2
  64..65      0.069    527.4    228.6   0.7310   0.0206   0.0281   10.8    6.4
  65..23      0.238    527.4    228.6   0.7310   0.0715   0.0978   37.7   22.4
  65..66      0.187    527.4    228.6   0.7310   0.0562   0.0769   29.6   17.6
  66..21      0.000    527.4    228.6   0.7310   0.0000   0.0000    0.0    0.0
  66..22      0.013    527.4    228.6   0.7310   0.0038   0.0052    2.0    1.2
  63..67      0.047    527.4    228.6   0.7310   0.0140   0.0192    7.4    4.4
  67..68      0.051    527.4    228.6   0.7310   0.0154   0.0211    8.1    4.8
  68..40      0.220    527.4    228.6   0.7310   0.0659   0.0902   34.8   20.6
  68..69      0.048    527.4    228.6   0.7310   0.0143   0.0195    7.5    4.5
  69..70      0.019    527.4    228.6   0.7310   0.0056   0.0076    2.9    1.7
  70..44      0.185    527.4    228.6   0.7310   0.0555   0.0759   29.3   17.3
  70..71      0.000    527.4    228.6   0.7310   0.0000   0.0000    0.0    0.0
  71..72      0.092    527.4    228.6   0.7310   0.0275   0.0377   14.5    8.6
  72..26      0.136    527.4    228.6   0.7310   0.0409   0.0559   21.6   12.8
  72..73      0.049    527.4    228.6   0.7310   0.0148   0.0202    7.8    4.6
  73..24      0.212    527.4    228.6   0.7310   0.0635   0.0868   33.5   19.8
  73..25      0.107    527.4    228.6   0.7310   0.0321   0.0439   16.9   10.0
  71..74      0.011    527.4    228.6   0.7310   0.0033   0.0045    1.7    1.0
  74..45      0.268    527.4    228.6   0.7310   0.0805   0.1101   42.4   25.2
  74..75      0.046    527.4    228.6   0.7310   0.0139   0.0191    7.3    4.4
  75..28      0.120    527.4    228.6   0.7310   0.0361   0.0494   19.1   11.3
  75..29      0.178    527.4    228.6   0.7310   0.0535   0.0732   28.2   16.7
  69..76      0.048    527.4    228.6   0.7310   0.0145   0.0198    7.6    4.5
  76..77      0.034    527.4    228.6   0.7310   0.0102   0.0140    5.4    3.2
  77..27      0.221    527.4    228.6   0.7310   0.0663   0.0907   35.0   20.7
  77..78      0.066    527.4    228.6   0.7310   0.0197   0.0269   10.4    6.1
  78..33      0.218    527.4    228.6   0.7310   0.0655   0.0896   34.5   20.5
  78..79      0.053    527.4    228.6   0.7310   0.0160   0.0219    8.4    5.0
  79..32      0.127    527.4    228.6   0.7310   0.0382   0.0523   20.2   12.0
  79..80      0.033    527.4    228.6   0.7310   0.0099   0.0135    5.2    3.1
  80..30      0.078    527.4    228.6   0.7310   0.0234   0.0320   12.3    7.3
  80..31      0.067    527.4    228.6   0.7310   0.0201   0.0275   10.6    6.3
  76..81      0.076    527.4    228.6   0.7310   0.0229   0.0313   12.1    7.1
  81..82      0.169    527.4    228.6   0.7310   0.0506   0.0692   26.7   15.8
  82..34      0.000    527.4    228.6   0.7310   0.0000   0.0000    0.0    0.0
  82..35      0.016    527.4    228.6   0.7310   0.0049   0.0067    2.6    1.5
  81..83      0.061    527.4    228.6   0.7310   0.0182   0.0249    9.6    5.7
  83..84      0.105    527.4    228.6   0.7310   0.0316   0.0432   16.6    9.9
  84..36      0.027    527.4    228.6   0.7310   0.0082   0.0112    4.3    2.6
  84..37      0.015    527.4    228.6   0.7310   0.0044   0.0060    2.3    1.4
  83..85      0.027    527.4    228.6   0.7310   0.0080   0.0109    4.2    2.5
  85..38      0.086    527.4    228.6   0.7310   0.0257   0.0351   13.5    8.0
  85..39      0.109    527.4    228.6   0.7310   0.0327   0.0447   17.2   10.2
  67..86      0.062    527.4    228.6   0.7310   0.0185   0.0253    9.8    5.8
  86..41      0.245    527.4    228.6   0.7310   0.0734   0.1005   38.7   23.0
  86..43      0.311    527.4    228.6   0.7310   0.0934   0.1278   49.3   29.2
  53..87      0.131    527.4    228.6   0.7310   0.0392   0.0537   20.7   12.3
  87..19      0.222    527.4    228.6   0.7310   0.0665   0.0910   35.1   20.8
  87..88      0.046    527.4    228.6   0.7310   0.0138   0.0189    7.3    4.3
  88..17      0.229    527.4    228.6   0.7310   0.0687   0.0940   36.2   21.5
  88..18      0.097    527.4    228.6   0.7310   0.0291   0.0398   15.3    9.1
上面的代码试图给树中的每一个结点加上codeml的计算结果.但是,程序设计的却有问题.%idlookup中有45个id索引.$tree->find_node(-id => $_)中的要求却大于45(因为45只代表树的叶结点,从上面的文件可知,包括内部结点,此树的结点共有88个),正如本程序所显示的一样"no node could be found for 53",找不到节点53.
因为这个程序主要是想知道处于正选择的位点,所以,干脆注释掉这一段程序,修改如下:

=mychange the number of branch_dnds and idlookup not match, so ....
        # These need to be added to the Node/branches

        while( my ($k,$v) = each %branch_dnds) {
         # we can probably do better by caching at some point

         my @nodes;
         for my $id ( split(/\.\./,$k ) ) {
            my @nodes_L = map { $tree->find_node(-id => $_) } @{$idlookup->{$id}};
            my $n = @nodes_L < 2 ? shift(@nodes_L) : $tree->get_lca(@nodes_L);
            if( ! $n ) {
             $self->warn("no node could be found for $id (no lca?)");
            }
            unless( $n->is_Leaf && $n->id) {
                $n->id($id);
                }
            push @nodes, $n;
         }
         my ($parent,$child) = @nodes;
         while ( my ($kk,$vv) = each %$v ) {
            $child->add_tag_value($kk,$vv);
         }
        }
=cut

运行结果如下:
-9241.335890
NEB
9 R 0.972 * 2.942
31 R 1.000 ** 2.998
32 C 0.936 2.869
33 N 1.000 ** 2.998
40 L 0.943 2.883
42 S 0.929 2.856
72 T 0.911 2.821
75 T 0.950 * 2.898
77 R 0.688 2.375
81 L 0.993 ** 2.983
101 P 0.995 ** 2.987
110 P 0.873 2.744
122 F 0.784 2.566
165 T 0.914 2.825
170 F 1.000 ** 2.997
182 V 1.000 ** 2.998
186 D 0.988 * 2.974
199 R 0.963 * 2.923
218 Y 0.879 2.755
219 D 0.970 * 2.938
221 T 0.999 ** 2.995
224 N 1.000 ** 2.998
225 R 1.000 ** 2.998
226 F 0.999 ** 2.997
232 H 0.608 2.214
252 C 0.863 2.725
BEB
9 R 0.965 * 2.974 0.616
31 R 1.000 ** 3.050 0.498
32 C 0.922 2.877 0.727
33 N 1.000 ** 3.051 0.498
40 L 0.932 2.903 0.703
42 S 0.919 2.874 0.734
72 T 0.890 2.803 0.790
75 T 0.925 2.875 0.718
77 R 0.605 2.125 0.983
81 L 0.992 ** 3.033 0.529
101 P 0.994 ** 3.039 0.521
110 P 0.812 2.606 0.896
122 F 0.742 2.467 0.966
165 T 0.836 2.654 0.864
170 F 1.000 ** 3.050 0.499
182 V 1.000 ** 3.050 0.498
186 D 0.985 * 3.016 0.555
199 R 0.940 2.910 0.683
218 Y 0.850 2.712   0.854
219 D 0.955 * 2.946 0.645
221 T 0.999 ** 3.047 0.503
224 N 1.000 ** 3.051 0.497
225 R 1.000 ** 3.051 0.497
226 F 0.999 ** 3.049 0.500
232 H 0.561 2.055 1.002
252 C 0.833 2.675 0.878
ps:附件是程序所需的文件
文件:ex4.tar.gz
大小:48KB
下载:下载
阅读(852) | 评论(0) | 转发(0) |
0

上一篇：[原创]又抓了bioperl中codeml.pm 模块的一个bug

下一篇：[ocrack 原创]bioperl 迁移到GitHub中了

给主人留下些什么吧！~~
关于我们 | 关于IT168 | 联系方式 | 广告合作 | 法律声明 | 免费注册

Copyright 2001-2010 ChinaUnix.net All Rights Reserved 北京皓辰网域网络信息技术有限公司. 版权所有

感谢所有关心和支持过ChinaUnix的朋友们

16024965号-6