← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to t/PAML.t

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
t/PAML.t
1
1
# This is -*-Perl-*- code
2
2
## Bioperl Test Harness Script for Modules
3
3
##
4
 
# $Id: PAML.t,v 1.12 2003/07/09 02:14:10 jason Exp $
 
4
# $Id: PAML.t,v 1.24.2.1 2006/11/08 17:25:55 sendu Exp $
5
5
 
6
6
# Before `make install' is performed this script should be runnable with
7
7
# `make test'. After `make install' it should work as `perl test.t'
9
9
use strict;
10
10
use vars qw($NUMTESTS $error);
11
11
 
12
 
 
13
12
BEGIN { 
14
13
    # to handle systems with no installed Test module
15
14
    # we include the t dir (where a copy of Test.pm is located)
21
20
    }
22
21
    use Test;
23
22
 
24
 
    $NUMTESTS = 116;
 
23
    $NUMTESTS = 192;
25
24
    plan tests => $NUMTESTS;
26
25
    eval { require IO::String; 
27
26
           require Bio::Tools::Phylo::PAML;}; 
28
 
    if( $@ ) { print STDERR "no IO string installed\n"; 
 
27
    if( $@ ) {
 
28
        print STDERR "no IO::String installed\n"; 
29
29
        $error = 1;
 
30
    }
 
31
}
 
32
 
 
33
END {
 
34
        foreach ( $Test::ntest .. $NUMTESTS ) {
 
35
                skip("Unable to run all of the PAML tests",1);
30
36
        }
31
37
}
32
38
 
33
 
END { 
34
 
    foreach ( $Test::ntest .. $NUMTESTS ) {
35
 
        skip("unable to run all of the PAML tests",1);
36
 
    }
37
 
}
38
 
 
39
39
 
40
40
exit(0) if( $error );
41
41
 
61
61
ok($NGmat->[0]->[1]->{'omega'}, 0.2507);
62
62
ok($NGmat->[0]->[1]->{'dN'}, 0.0863);
63
63
ok($NGmat->[0]->[1]->{'dS'}, 0.3443);
 
64
 
64
65
ok($NGmat->[2]->[3]->{'omega'}, 0.2178);
65
66
ok($NGmat->[2]->[3]->{'dN'}, 0.1348);
66
67
ok($NGmat->[2]->[3]->{'dS'}, 0.6187);
67
68
 
68
 
ok($MLmat->[0]->[1]->{'omega'}, 0.1948);
 
69
ok($MLmat->[0]->[1]->{'omega'}, 0.19479);
69
70
ok($MLmat->[0]->[1]->{'dN'}, 0.0839);
70
71
ok($MLmat->[0]->[1]->{'dS'}, 0.4309);
71
72
ok($MLmat->[0]->[1]->{'lnL'}, -1508.607268);
72
 
ok($MLmat->[2]->[3]->{'omega'}, 0.1611);
 
73
ok($MLmat->[0]->[1]->{'t'}, 0.47825);
 
74
ok($MLmat->[0]->[1]->{'kappa'}, 2.29137);
 
75
 
 
76
ok($MLmat->[2]->[3]->{'omega'}, 0.16114);
73
77
ok($MLmat->[2]->[3]->{'dN'}, 0.1306);
74
78
ok($MLmat->[2]->[3]->{'dS'}, 0.8105);
75
79
ok($MLmat->[2]->[3]->{'lnL'},-1666.440696);
 
80
ok($MLmat->[2]->[3]->{'t'}, 0.85281);
 
81
ok($MLmat->[2]->[3]->{'kappa'}, 2.21652);
76
82
 
77
83
my @codonposfreq = $result->get_codon_pos_basefreq();
78
84
ok($codonposfreq[0]->{'A'}, 0.23579);
244
250
ok($firstsite->[1], 'L');
245
251
ok($firstsite->[2], 0.6588);
246
252
 
 
253
# codeml NSSites parsing
 
254
# for M0 model
 
255
 
 
256
my $codeml_m0 = new Bio::Tools::Phylo::PAML
 
257
    (-file => Bio::Root::IO->catfile(qw/t data M0.mlc/));
 
258
ok($codeml_m0);
 
259
my $result_m0 = $codeml_m0->next_result;
 
260
my ($nssite_m0,$nssite_m1) = $result_m0->get_NSSite_results;
 
261
ok($nssite_m0->num_site_classes,1);
 
262
my $class_m0 = $nssite_m0->dnds_site_classes;
 
263
ok($class_m0->{q/p/}->[0],q/1.00000/);
 
264
ok($class_m0->{q/w/}->[0],0.09213);
 
265
 
 
266
ok($nssite_m0->model_num, "0");
 
267
@trees= $nssite_m0->get_trees;
 
268
ok (@trees , 1 ); 
 
269
# model 0
 
270
ok($trees[0]->score, -30.819156);
 
271
ok($nssite_m1->model_num, "1");
 
272
@trees= $nssite_m1->get_trees;
 
273
ok($trees[0]->score, -30.819157);
 
274
 
 
275
# test BASEML
 
276
# pairwise first
 
277
 
 
278
my $baseml_p = Bio::Tools::Phylo::PAML->new
 
279
    (-file => Bio::Root::IO->catfile(qw(t data baseml.pairwise)));
 
280
ok($baseml_p);
 
281
my $baseml = $baseml_p->next_result;
 
282
my @b_seqs =  $baseml->get_seqs;
 
283
ok($b_seqs[0]->seq, 'GTAGAGTACTTT');
 
284
ok($b_seqs[1]->seq, 'GTAAGAGACGAT');
 
285
 
 
286
my @otus = map { $_->display_id } @b_seqs;
 
287
ok(scalar @otus, 3);
 
288
my $ntfreq = $baseml->get_NTFreqs;
 
289
ok($ntfreq);
 
290
ok($ntfreq->{$otus[0]}->{'A'}, '0.3333');
 
291
ok($ntfreq->{$otus[1]}->{'G'}, '0.2105');
 
292
my $kappaM = $baseml->get_KappaMatrix;
 
293
ok($kappaM);
 
294
ok($kappaM->get_entry($otus[1],$otus[0]), '0.3240');
 
295
ok($kappaM->get_entry($otus[0],$otus[1]), 
 
296
   $kappaM->get_entry($otus[1],$otus[0]));
 
297
ok($kappaM->get_entry($otus[1],$otus[2]), '0.1343');
 
298
my $alphaM = $baseml->get_AlphaMatrix;
 
299
ok($alphaM);
 
300
ok($alphaM->get_entry($otus[1],$otus[0]), '9.3595');
 
301
ok($alphaM->get_entry($otus[0],$otus[1]), 
 
302
   $alphaM->get_entry($otus[1],$otus[0]));
 
303
ok($alphaM->get_entry($otus[1],$otus[2]), '1.1101');
 
304
ok($alphaM->get_entry($otus[0],$otus[2]), '33.1197');
 
305
 
 
306
# codeml NSSites parsing
 
307
# for only 1 model
 
308
 
 
309
my $codeml_single = new Bio::Tools::Phylo::PAML
 
310
    (-file => Bio::Root::IO->catfile(qw/t data singleNSsite.mlc/));
 
311
ok($codeml_single);
 
312
my $result_single = $codeml_single->next_result;
 
313
my ($nssite_single) = $result_single->get_NSSite_results;
 
314
ok($nssite_single->num_site_classes,q/3/);
 
315
ok($nssite_single->kappa, q/5.28487/);
 
316
ok($nssite_single->likelihood,q/-30.819156/);
 
317
 
 
318
ok($baseml->get_stat('loglikelihood'),-110.532715);
 
319
ok($baseml->get_stat('constant_sites'),46);
 
320
ok($baseml->get_stat('constant_sites_percentage'),'80.70');
 
321
ok($baseml->model,'HKY85 dGamma (ncatG=5)');
 
322
 
 
323
# user trees
 
324
$baseml_p = Bio::Tools::Phylo::PAML->new
 
325
    (-file => Bio::Root::IO->catfile(qw(t data baseml.usertree)));
 
326
$baseml = $baseml_p->next_result;
 
327
 
 
328
@trees = $baseml->get_trees;
 
329
ok(@trees, 1);
 
330
ok($trees[0]->score, -129.328757);
 
331
 
 
332
# codeml NSSites parsing
 
333
# for branch site model/clade model
 
334
 
 
335
my $codeml_bs = new Bio::Tools::Phylo::PAML
 
336
    (-file => Bio::Root::IO->catfile(qw/t data branchSite.mlc/));
 
337
ok($codeml_bs);
 
338
my $result_bs = $codeml_bs->next_result;
 
339
my ($nssite_bs) = $result_bs->get_NSSite_results;
 
340
ok($nssite_bs->num_site_classes,q/4/);
 
341
my $class_bs = $nssite_bs->dnds_site_classes;
 
342
ok($class_bs->{q/p/}->[1],q/0.65968/);
 
343
ok($class_bs->{q/w/}->[1]->{q/background/},q/0.00000/);
 
344
ok($class_bs->{q/w/}->[2]->{q/foreground/},q/999.00000/);
 
345
 
 
346
# Let's parse the RST file
 
347
 
 
348
my $paml = Bio::Tools::Phylo::PAML->new
 
349
    (-file => Bio::Root::IO->catfile(qw(t data codeml_lysozyme mlc)),
 
350
     -dir  => Bio::Root::IO->catfile(qw(t data codeml_lysozyme)));
 
351
 
 
352
$result = $paml->next_result;
 
353
 
 
354
my ($rst) = grep {$_->id eq 'node#8'} $result->get_rst_seqs;
 
355
ok($rst);
 
356
ok($rst->seq, join('',qw(
 
357
AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACTCTGAAAAGATTGGGACTGGATGGCTAC
 
358
AGGGGAATCAGCCTAGCAAACTGGATGTGTTTGGCCAAATGGGAGAGTGATTATAACACA
 
359
CGAGCTACAAACTACAATCCTGGAGACCAAAGCACTGATTATGGGATATTTCAGATCAAT
 
360
AGCCACTACTGGTGTAATAATGGCAAAACCCCAGGAGCAGTTAATGCCTGTCATATATCC
 
361
TGCAATGCTTTGCTGCAAGATAACATCGCTGATGCTGTAGCTTGTGCAAAGAGGGTTGTC
 
362
CGTGATCCACAAGGCATTAGAGCATGGGTGGCATGGAGAAATCATTGTCAAAACAGAGAT
 
363
GTCAGTCAGTATGTTCAAGGTTGTGGAGTG)),
 
364
   'node#8 reconstructed seq');
 
365
 
 
366
my ($first_tree) = $result->get_rst_trees;
 
367
my ($node) = $first_tree->find_node(-id => '5_Mmu_rhesus');
 
368
my @changes = $node->get_tag_values('changes');
 
369
my ($site) = grep { $_->{'site'} == 94 } @changes;
 
370
ok($site->{'anc_aa'}, 'A','ancestral AA');
 
371
ok($site->{'anc_prob'}, '0.947','ancestral AA prob');
 
372
ok($site->{'derived_aa'}, 'T','derived AA');
 
373
 
 
374
($node) = $first_tree->find_node(-id => '12');
 
375
@changes = $node->get_tag_values('changes');
 
376
($site) = grep { $_->{'site'} == 88 } @changes;
 
377
ok($site->{'anc_aa'}, 'N','ancestral AA');
 
378
ok($site->{'anc_prob'}, '0.993','ancestral AA prob');
 
379
ok($site->{'derived_aa'}, 'D','derived AA');
 
380
ok($site->{'derived_prob'}, '0.998','derived AA prob');
 
381
 
 
382
my $persite = $result->get_rst_persite;
 
383
# minus 1 because we have shifted so that array index matches site number
 
384
# there are 130 sites in this seq file
 
385
ok(scalar @$persite -1, $result->patterns->{'-ls'}); 
 
386
# let's score site 1
 
387
$site = $persite->[2]; 
 
388
# so site 2, node 2 (extant)
 
389
ok($site->[2]->{'codon'}, 'GTC');
 
390
ok($site->[2]->{'aa'}, 'V');
 
391
# site 2, node 3
 
392
ok($site->[3]->{'codon'}, 'ATC');
 
393
ok($site->[3]->{'aa'}, 'I');
 
394
 
 
395
# ancestral node 9
 
396
ok($site->[9]->{'codon'}, 'GTC');
 
397
ok($site->[9]->{'aa'},    'V');
 
398
ok($site->[9]->{'prob'},  '1.000');
 
399
ok($site->[9]->{'Yang95_aa'},'V');
 
400
ok($site->[9]->{'Yang95_aa_prob'},'1.000');
 
401
 
 
402
# ancestral node 10
 
403
ok($site->[10]->{'codon'}, 'ATC');
 
404
ok($site->[10]->{'aa'},    'I');
 
405
ok($site->[10]->{'prob'},  '0.992');
 
406
ok($site->[10]->{'Yang95_aa'},'I');
 
407
ok($site->[10]->{'Yang95_aa_prob'},'0.992');
 
408
 
 
409
 
 
410
## PAML 3.15
 
411
$paml = Bio::Tools::Phylo::PAML->new(-file => Bio::Root::IO->catfile(qw(t data codeml315.mlc)) );
 
412
$result = $paml->next_result;
 
413
 
 
414
ok($result->model, 'One dN/dS ratio');
 
415
ok($result->version, qr'3\.15');
 
416
$MLmat = $result->get_MLmatrix;
 
417
$NGmat = $result->get_NGmatrix;
 
418
 
 
419
ok($NGmat->[0]->[1]->{'omega'}, 0.2264);
 
420
ok($NGmat->[0]->[1]->{'dN'}, 0.0186);
 
421
ok($NGmat->[0]->[1]->{'dS'}, 0.0821);
 
422
 
 
423
ok($MLmat->[0]->[1]->{'omega'}, 0.32693);
 
424
ok($MLmat->[0]->[1]->{'dN'}, '0.0210');
 
425
ok($MLmat->[0]->[1]->{'dS'}, 0.0644);