← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/bioperl/oneiric

« back to all changes in this revision

Viewing changes to t/PAML.t

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
t/PAML.t
 
1
# This is -*-Perl-*- code
 
2
## Bioperl Test Harness Script for Modules
 
3
##
 
4
# $Id: PAML.t,v 1.12 2003/07/09 02:14:10 jason Exp $
 
5
 
 
6
# Before `make install' is performed this script should be runnable with
 
7
# `make test'. After `make install' it should work as `perl test.t'
 
8
 
 
9
use strict;
 
10
use vars qw($NUMTESTS $error);
 
11
 
 
12
 
 
13
BEGIN { 
 
14
    # to handle systems with no installed Test module
 
15
    # we include the t dir (where a copy of Test.pm is located)
 
16
    # as a fallback
 
17
    eval { require Test; };
 
18
    $error = 0;
 
19
    if( $@ ) {
 
20
        use lib 't';
 
21
    }
 
22
    use Test;
 
23
 
 
24
    $NUMTESTS = 116;
 
25
    plan tests => $NUMTESTS;
 
26
    eval { require IO::String; 
 
27
           require Bio::Tools::Phylo::PAML;}; 
 
28
    if( $@ ) { print STDERR "no IO string installed\n"; 
 
29
        $error = 1;
 
30
        }
 
31
}
 
32
 
 
33
END { 
 
34
    foreach ( $Test::ntest .. $NUMTESTS ) {
 
35
        skip("unable to run all of the PAML tests",1);
 
36
    }
 
37
}
 
38
 
 
39
 
 
40
exit(0) if( $error );
 
41
 
 
42
my $testnum;
 
43
my $verbose = 0;
 
44
 
 
45
## End of black magic.
 
46
##
 
47
## Insert additional test code below but remember to change
 
48
## the print "1..x\n" in the BEGIN block to reflect the
 
49
## total number of tests that will be run. 
 
50
use Bio::Root::IO;
 
51
 
 
52
my $inpaml = new Bio::Tools::Phylo::PAML(-file => Bio::Root::IO->catfile(qw(t data codeml.mlc)));
 
53
ok($inpaml);
 
54
my $result = $inpaml->next_result;
 
55
ok($result);
 
56
ok($result->model, 'several dN/dS ratios for branches');
 
57
ok($result->version, qr'3\.12');
 
58
my $MLmat = $result->get_MLmatrix;
 
59
my $NGmat = $result->get_NGmatrix;
 
60
 
 
61
ok($NGmat->[0]->[1]->{'omega'}, 0.2507);
 
62
ok($NGmat->[0]->[1]->{'dN'}, 0.0863);
 
63
ok($NGmat->[0]->[1]->{'dS'}, 0.3443);
 
64
ok($NGmat->[2]->[3]->{'omega'}, 0.2178);
 
65
ok($NGmat->[2]->[3]->{'dN'}, 0.1348);
 
66
ok($NGmat->[2]->[3]->{'dS'}, 0.6187);
 
67
 
 
68
ok($MLmat->[0]->[1]->{'omega'}, 0.1948);
 
69
ok($MLmat->[0]->[1]->{'dN'}, 0.0839);
 
70
ok($MLmat->[0]->[1]->{'dS'}, 0.4309);
 
71
ok($MLmat->[0]->[1]->{'lnL'}, -1508.607268);
 
72
ok($MLmat->[2]->[3]->{'omega'}, 0.1611);
 
73
ok($MLmat->[2]->[3]->{'dN'}, 0.1306);
 
74
ok($MLmat->[2]->[3]->{'dS'}, 0.8105);
 
75
ok($MLmat->[2]->[3]->{'lnL'},-1666.440696);
 
76
 
 
77
my @codonposfreq = $result->get_codon_pos_basefreq();
 
78
ok($codonposfreq[0]->{'A'}, 0.23579);
 
79
ok($codonposfreq[0]->{'T'}, 0.14737);
 
80
ok($codonposfreq[1]->{'C'}, 0.25123);
 
81
ok($codonposfreq[2]->{'G'}, 0.32842);
 
82
 
 
83
# AAML parsing - Empirical model
 
84
$inpaml = new Bio::Tools::Phylo::PAML(-file => Bio::Root::IO->catfile
 
85
                                      (qw(t data aaml.mlc)));
 
86
 
 
87
ok($inpaml);
 
88
$result = $inpaml->next_result;
 
89
ok($result);
 
90
ok($result->model, 'Empirical (wag.dat)');
 
91
my @trees = $result->get_trees;
 
92
ok(@trees, 1);
 
93
ok($trees[0]->score, -1042.768973);
 
94
 
 
95
ok((scalar grep { $_->is_Leaf } $trees[0]->get_nodes), $result->get_seqs);
 
96
 
 
97
my $aadistmat = $result->get_AADistMatrix();
 
98
ok($aadistmat);
 
99
ok($aadistmat->get_entry('Cow', 'Horse'), 0.5462);
 
100
ok($aadistmat->get_entry('Baboon', 'Langur'), 0.1077);
 
101
 
 
102
my %aafreq = %{$result->get_AAFreqs()};
 
103
ok(%aafreq);
 
104
ok($aafreq{'Human'}->{'N'}, 0.0769);
 
105
ok($aafreq{'Human'}->{'R'}, 0.1077);
 
106
 
 
107
my %ratfreqs = %{$result->get_AAFreqs('Rat')};
 
108
ok($ratfreqs{'R'},0.0923);
 
109
ok($ratfreqs{'F'},0.0154);
 
110
my %avgfreqs = %{$result->get_AAFreqs('Average')};
 
111
ok($avgfreqs{'Q'},0.0411);
 
112
 
 
113
ok($result->get_AAFreqs('Average')->{'I'},0.0424);
 
114
 
 
115
my $patterns = $result->patterns;
 
116
my @pat = @{$patterns->{'-patterns'}};
 
117
ok(scalar @pat, 98);
 
118
ok($patterns->{'-ns'}, 6);
 
119
ok($patterns->{'-ls'}, 130);
 
120
 
 
121
ok((sort $result->get_stat_names)[0], 'constant_sites');
 
122
ok($result->get_stat('constant_sites'), 46);
 
123
ok($result->get_stat('constant_sites_percentage'), 35.38);
 
124
 
 
125
# AAML parsing - pairwise model
 
126
$inpaml = new Bio::Tools::Phylo::PAML(-file => Bio::Root::IO->catfile
 
127
                                      (qw(t data aaml_pairwise.mlc)));
 
128
 
 
129
ok($inpaml);
 
130
$result = $inpaml->next_result;
 
131
ok($result);
 
132
ok($result->model, 'Empirical_F (wag.dat)');
 
133
ok($result->get_stat('loglikelihood'),-1189.106658);
 
134
ok($result->get_stat('constant_sites'), 170);
 
135
ok($result->get_stat('constant_sites_percentage'), 59.65);
 
136
 
 
137
ok($result->get_AAFreqs('Average')->{'R'},0.0211);
 
138
ok($result->get_AAFreqs('rabbit')->{'L'},0.1228);
 
139
 
 
140
$aadistmat = $result->get_AADistMatrix();
 
141
ok($aadistmat);
 
142
ok($aadistmat->get_entry('rabbit', 'marsupial'), 0.2877);
 
143
ok($aadistmat->get_entry('human', 'goat-cow'), 0.1439);
 
144
 
 
145
$aadistmat = $result->get_AAMLDistMatrix();
 
146
ok($aadistmat);
 
147
ok($aadistmat->get_entry('rabbit', 'marsupial'), 0.3392);
 
148
ok($aadistmat->get_entry('human', 'goat-cow'), 0.1551);
 
149
 
 
150
my @seqs = $result->get_seqs;
 
151
ok($seqs[0]->display_id, 'human');
 
152
 
 
153
# YN00 parsing, pairwise Ka/Ks from Yang & Nielsen 2000
 
154
$inpaml = new Bio::Tools::Phylo::PAML(-file => Bio::Root::IO->catfile
 
155
                                      (qw(t data yn00.mlc)));
 
156
 
 
157
ok($inpaml);
 
158
$result = $inpaml->next_result;
 
159
 
 
160
ok($result);
 
161
$MLmat = $result->get_MLmatrix;
 
162
$NGmat = $result->get_NGmatrix;
 
163
 
 
164
ok($NGmat->[0]->[1]->{'omega'}, 0.251);
 
165
ok($NGmat->[0]->[1]->{'dN'}, 0.0863);
 
166
ok($NGmat->[0]->[1]->{'dS'}, 0.3443);
 
167
ok($NGmat->[2]->[3]->{'omega'}, 0.218);
 
168
ok($NGmat->[2]->[3]->{'dN'}, 0.1348);
 
169
ok($NGmat->[2]->[3]->{'dS'}, 0.6187);
 
170
 
 
171
ok($MLmat->[0]->[1]->{'omega'}, 0.1625);
 
172
ok($MLmat->[0]->[1]->{'dN'}, 0.0818);
 
173
ok($MLmat->[0]->[1]->{'dS'}, 0.5031);
 
174
ok($MLmat->[2]->[3]->{'omega'}, 0.1262);
 
175
ok($MLmat->[2]->[3]->{'dN'}, 0.1298);
 
176
ok($MLmat->[2]->[3]->{'dN_SE'}, 0.0149);
 
177
ok($MLmat->[2]->[3]->{'dS'}, 1.0286);
 
178
ok($MLmat->[2]->[3]->{'dS_SE'}, 0.2614);
 
179
 
 
180
# codeml NSSites parsing
 
181
 
 
182
$inpaml = new Bio::Tools::Phylo::PAML
 
183
    (-file => Bio::Root::IO->catfile(qw(t data codeml_nssites.mlc)));
 
184
 
 
185
ok($inpaml);
 
186
$result = $inpaml->next_result;
 
187
 
 
188
ok($result);
 
189
ok($result->model, 'One dN/dS ratio dGamma (ncatG=11)');
 
190
ok($result->version, 'paml 3.13, August 2002');
 
191
$NGmat = $result->get_NGmatrix;
 
192
ok($NGmat);
 
193
 
 
194
ok($NGmat->[0]->[1]->{'omega'}, 0.2782);
 
195
ok($NGmat->[0]->[1]->{'dN'}, 0.0133);
 
196
ok($NGmat->[0]->[1]->{'dS'}, 0.0478);
 
197
ok($NGmat->[1]->[2]->{'omega'}, 1.1055);
 
198
ok($NGmat->[1]->[2]->{'dN'}, 0.0742);
 
199
ok($NGmat->[1]->[2]->{'dS'}, 0.0671);
 
200
          # this is
 
201
          #   model num  description
 
202
          #   kappa   log-likelihood tree length time used
 
203
          #   shape   alpha/gamma r          f
 
204
my @tstr = ([qw(0 one-ratio 0
 
205
                4.54006 -906.017440    0.55764
 
206
                )],
 
207
            [qw(1 neutral 2
 
208
                4.29790 -902.503869    0.56529
 
209
                )],
 
210
            [qw(2 selection 3 
 
211
                5.12250 -900.076500    0.6032
 
212
                )],
 
213
             );
 
214
my $iter = 0;
 
215
my $lastmodel;
 
216
foreach my $model ( $result->get_NSSite_results ) {    
 
217
    my $i = 0;
 
218
    my $r = shift @tstr;
 
219
    ok($model->model_num, $r->[$i++]);
 
220
    ok($model->model_description, qr/$r->[$i++]/);
 
221
    ok($model->num_site_classes,$r->[$i++]);
 
222
    my $tree = $model->next_tree;
 
223
    ok($model->kappa, $r->[$i++]);
 
224
    ok($model->likelihood,$r->[$i]);
 
225
    ok($tree->score, $r->[$i++]);
 
226
    ok($tree->total_branch_length, $r->[$i++]);
 
227
    if( $iter == 0 ) {
 
228
        my $params = $model->shape_params;
 
229
        ok($params->{'shape'}, 'alpha');
 
230
        ok($params->{'gamma'},   '0.50000');
 
231
        ok($params->{'r'}->[0], '1.00000');
 
232
        ok($params->{'f'}->[0], '1.00000');
 
233
    } elsif( $iter == 2 ) {
 
234
        my $class = $model->dnds_site_classes;
 
235
        ok($class->{'p'}->[0], '0.38160');
 
236
        ok($class->{'w'}->[1], '1.00000');
 
237
    }
 
238
    $iter++;
 
239
    $lastmodel = $model;
 
240
}
 
241
 
 
242
my ($firstsite) = $lastmodel->get_pos_selected_sites;
 
243
ok($firstsite->[0], 15);
 
244
ok($firstsite->[1], 'L');
 
245
ok($firstsite->[2], 0.6588);
 
246