t/AlignStats.t

   1 # -*-Perl-*-
   2 ## Bioperl Test Harness Script for Modules
   3 ## $Id$
   4
   5 # Before `make install' is performed this script should be runnable with
   6 # `make test'. After `make install' it should work as `perl test.t'
   7
   8 my $error = 0;
   9
  10 use strict;
  11 BEGIN {
  12     # to handle systems with no installed Test module
  13     # we include the t dir (where a copy of Test.pm is located)
  14     # as a fallback
  15     eval { require Test; };
  16     if( $@ ) {
  17         use lib 't';
  18     }
  19
  20     use Test;
  21     plan tests => 36;
  22 }
  23
  24 if( $error == 1 ) {
  25     exit(0);
  26 }
  27
  28 my $debug = -1;
  29
  30 use Bio::Align::DNAStatistics;
  31 use Bio::Align::ProteinStatistics;
  32 use Bio::AlignIO;
  33 use Bio::Root::IO;
  34
  35 my $in = new Bio::AlignIO(-format => 'emboss',
  36                           -file   => Bio::Root::IO->catfile('t', 'data',
  37                                                             'insulin.water'));
  38 my $aln = $in->next_aln();
  39 ok($aln);
  40 my $stats = new Bio::Align::DNAStatistics(-verbose => $debug);
  41 ok( $stats->transversions($aln),4);
  42 ok( $stats->transitions($aln),9);
  43 ok( $stats->pairwise_stats->number_of_gaps($aln),21);
  44 ok( $stats->pairwise_stats->number_of_comparable_bases($aln),173);
  45 ok( $stats->pairwise_stats->number_of_differences($aln),13);
  46
  47 my $d = $stats->distance(-align => $aln,
  48                          -method=> 'f81');
  49 ok(  $d->get_entry('hs_insulin','seq2'), '0.07918');
  50
  51 $d = $stats->distance(-align=> $aln,
  52                       -method => 'JC');
  53 ok( $d->get_entry('hs_insulin','seq2'), '0.07918');
  54
  55 $d = $stats->distance(-align=> $aln,
  56                       -method => 'Kimura');
  57 ok( $d->get_entry('hs_insulin','seq2'), '0.07984');
  58
  59 $d = $stats->distance(-align=> $aln,
  60                       -method => 'TajimaNei');
  61 ok( $d->get_entry('seq2','hs_insulin'), '0.08106');
  62
  63 $d = $stats->distance(-align=> $aln,
  64                       -method => 'Tamura');
  65 ok( $d->get_entry('seq2','hs_insulin'), '0.08037');
  66
  67 #$d =  $stats->distance(-align => $aln,
  68 #                      -method => 'JinNei');
  69 #ok( $d->get_entry('seq2','hs_insulin'), 0.0850);
  70
  71 $in = new Bio::AlignIO(-format => 'clustalw',
  72                        -file   => Bio::Root::IO->catfile('t','data',
  73                                                          'hs_owlmonkey.aln'));
  74
  75 $aln = $in->next_aln();
  76 ok($aln);
  77
  78 ok( $stats->transversions($aln),10);
  79 ok( $stats->transitions($aln),17);
  80 ok( $stats->pairwise_stats->number_of_gaps($aln),19);
  81 ok( $stats->pairwise_stats->number_of_comparable_bases($aln),170);
  82 ok( $stats->pairwise_stats->number_of_differences($aln),27);
  83
  84 # now test the distance calculations
  85 $d = $stats->distance(-align => $aln, -method => 'jc');
  86 ok( $d->get_entry('human','owlmonkey'), 0.17847);
  87
  88 $d = $stats->distance(-align => $aln,
  89                       -method=> 'f81');
  90 ok(  $d->get_entry('human','owlmonkey'), '0.17847');
  91
  92 $d = $stats->distance(-align => $aln, -method => 'uncorrected');
  93 ok( $d->get_entry('human','owlmonkey'), 0.15882);
  94
  95 $d =  $stats->distance(-align => $aln, -method => 'Kimura');
  96 ok( $d->get_entry('human','owlmonkey'), 0.18105);
  97
  98 $d =  $stats->distance(-align => $aln, -method => 'TajimaNei');
  99 ok( $d->get_entry('human','owlmonkey'), 0.18489);
 100
 101 $d =  $stats->distance(-align => $aln,
 102                        -method => 'Tamura');
 103
 104 ok( $d->get_entry('human','owlmonkey'), 0.18333);
 105 #$d =  $stats->distance(-align => $aln,
 106 #                      -method => 'JinNei');
 107 #ok( $d->get_entry('human','owlmonkey'), 0.2079);
 108
 109 ### now test Nei_gojobori methods ##
 110 $in = Bio::AlignIO->new(-format => 'fasta',
 111                         -file   => Bio::Root::IO->catfile('t','data',
 112                                                           'nei_gojobori_test.aln'));
 113 my $alnobj = $in->next_aln();
 114 ok($alnobj);
 115 my $result = $stats->calc_KaKs_pair($alnobj, 'seq1', 'seq2');
 116 ok (sprintf ("%.1f", $result->[0]{'S'}), 40.5);
 117 ok (sprintf ("%.1f", $result->[0]{'z_score'}), '4.5');
 118 $result = $stats->calc_all_KaKs_pairs($alnobj);
 119 ok (int( $result->[1]{'S'}), 41);
 120 ok (int( $result->[1]{'z_score'}), 4);
 121 $result = $stats->calc_average_KaKs($alnobj, 100);
 122 ok (sprintf ("%.4f", $result->{'D_n'}), 0.1628);
 123
 124
 125 # now test Protein Distances
 126 my $pstats = Bio::Align::ProteinStatistics->new();
 127 $in = Bio::AlignIO->new(-format => 'clustalw',
 128                         -file   => Bio::Root::IO->catfile('t','data',
 129                                                           'testaln.aln'));
 130 $alnobj = $in->next_aln();
 131 ok($alnobj);
 132 $result = $pstats->distance(-method => 'Kimura',
 133                             -align  => $alnobj);
 134 ok($result);
 135
 136 ok ($result->get_entry('P84139','P814153'),   '0.01443');
 137 ok ($result->get_entry('P841414','P851414'),  '0.01686');
 138 ok ($result->get_entry('P84139','P851414'),   '3.58352');
 139
 140 my $seq = Bio::Seq->new(-id=>'NOT3MUL', -seq=>'gatac');
 141 ok($seq);
 142 eval {
 143   Bio::Align::DNAStatistics->count_syn_sites($seq);
 144 };
 145 ok($@ =~ m/not integral number of codons/);
 146
 147