t/Align/AlignStats.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use Bio::Root::Test;
   8
   9     test_begin(-tests => 45);
  10
  11         use_ok('Bio::Align::DNAStatistics');
  12         use_ok('Bio::Align::ProteinStatistics');
  13         use_ok('Bio::AlignIO');
  14 }
  15
  16 my $debug = test_debug();
  17
  18 my $in = Bio::AlignIO->new(-format => 'emboss',
  19                           -file   => test_input_file('insulin.water'));
  20 my $aln = $in->next_aln();
  21 isa_ok($aln, 'Bio::Align::AlignI');
  22 my $stats = Bio::Align::DNAStatistics->new(-verbose => $debug);
  23 is( $stats->transversions($aln),4);
  24 is( $stats->transitions($aln),9);
  25 is( $stats->pairwise_stats->number_of_gaps($aln),21);
  26 is( $stats->pairwise_stats->number_of_comparable_bases($aln),173);
  27 is( $stats->pairwise_stats->number_of_differences($aln),13);
  28 is( $stats->pairwise_stats->score_nuc($aln), 224);
  29 is( $stats->pairwise_stats->score_nuc( -aln => $aln, -match => 1,
  30   -mismatch => -1, -gap_open => -1, -gap_ext => -1), 126);
  31
  32 my $d = $stats->distance(-align => $aln,
  33                          -method=> 'f81');
  34 is(  $d->get_entry('hs_insulin','seq2'), '0.07918');
  35
  36 $d = $stats->distance(-align=> $aln,
  37                       -method => 'JC');
  38 is( $d->get_entry('hs_insulin','seq2'), '0.07918');
  39
  40 $d = $stats->distance(-align=> $aln,
  41                       -method => 'Kimura');
  42 is( $d->get_entry('hs_insulin','seq2'), '0.07984');
  43
  44 $d = $stats->distance(-align=> $aln,
  45                       -method => 'TajimaNei');
  46 is( $d->get_entry('seq2','hs_insulin'), '0.08106');
  47
  48 $d = $stats->distance(-align=> $aln,
  49                       -method => 'Tamura');
  50 is( $d->get_entry('seq2','hs_insulin'), '0.08037');
  51
  52 #$d =  $stats->distance(-align => $aln,
  53 #                      -method => 'JinNei');
  54 #is( $d->get_entry('seq2','hs_insulin'), 0.0850);
  55
  56 $in = Bio::AlignIO->new(-format => 'clustalw',
  57                        -file   => test_input_file('hs_owlmonkey.aln'));
  58
  59 $aln = $in->next_aln();
  60 isa_ok($aln,'Bio::Align::AlignI');
  61
  62 is( $stats->transversions($aln),10);
  63 is( $stats->transitions($aln),17);
  64 is( $stats->pairwise_stats->number_of_gaps($aln),19);
  65 is( $stats->pairwise_stats->number_of_comparable_bases($aln),170);
  66 is( $stats->pairwise_stats->number_of_differences($aln),27);
  67 is( $stats->pairwise_stats->score_nuc($aln), 134);
  68 is( $stats->pairwise_stats->score_nuc( -aln => $aln, -match => 1,
  69   -mismatch => -1, -gap_open => -1, -gap_ext => -1), 97);
  70
  71 # now test the distance calculations
  72 $d = $stats->distance(-align => $aln, -method => 'jc');
  73 is( $d->get_entry('human','owlmonkey'), 0.17847);
  74
  75 $d = $stats->distance(-align => $aln,
  76                           -method=> 'f81');
  77 is(  $d->get_entry('human','owlmonkey'), '0.17847');
  78
  79 $d = $stats->distance(-align => $aln, -method => 'uncorrected');
  80 is( $d->get_entry('human','owlmonkey'), 0.15882);
  81
  82 $d =  $stats->distance(-align => $aln, -method => 'Kimura');
  83 is( $d->get_entry('human','owlmonkey'), 0.18105);
  84
  85 $d =  $stats->distance(-align => $aln, -method => 'TajimaNei');
  86 is( $d->get_entry('human','owlmonkey'), 0.18489);
  87
  88 $d =  $stats->distance(-align => $aln,
  89                            -method => 'Tamura');
  90
  91 is( $d->get_entry('human','owlmonkey'), 0.18333);
  92 #$d =  $stats->distance(-align => $aln,
  93 #                      -method => 'JinNei');
  94 #is( $d->get_entry('human','owlmonkey'), 0.2079);
  95
  96 ### now test Nei_gojobori methods, hiding the expected warnings so we can
  97 # avoid printing them ###
  98 $stats->verbose($debug ? $debug : -1);
  99 my ($alnobj, $result);
 100 $in = Bio::AlignIO->new(-format => 'fasta',
 101                         -file   => test_input_file('nei_gojobori_test.aln'));
 102 $alnobj = $in->next_aln();
 103 isa_ok($alnobj,'Bio::Align::AlignI');
 104 $result = $stats->calc_KaKs_pair($alnobj, 'seq1', 'seq2');
 105 is (sprintf ("%.1f", $result->[0]{'S'}), 40.5);
 106 is (sprintf ("%.1f", $result->[0]{'z_score'}), '4.5');
 107 $result = $stats->calc_all_KaKs_pairs($alnobj);
 108 is (int( $result->[1]{'S'}), 41);
 109 is (int( $result->[1]{'z_score'}), 4);
 110 $result = $stats->calc_average_KaKs($alnobj, 100);
 111 is (sprintf ("%.4f", $result->{'D_n'}), 0.1628);
 112 $stats->verbose($debug);
 113
 114 # now test Protein Distances
 115 my $pstats = Bio::Align::ProteinStatistics->new();
 116 $in = Bio::AlignIO->new(-format => 'clustalw',
 117                         -file   => test_input_file('testaln.clustalw'));
 118 $alnobj = $in->next_aln();
 119 isa_ok($alnobj,'Bio::Align::AlignI');
 120 $result = $pstats->distance(-method => 'Kimura',
 121                             -align  => $alnobj);
 122 isa_ok($result, 'Bio::Matrix::PhylipDist');
 123
 124 is ($result->get_entry('P84139','P814153'),   '0.01443');
 125 is ($result->get_entry('P841414','P851414'),  '0.01686');
 126 is ($result->get_entry('P84139','P851414'),   '3.58352');
 127
 128 my $seq = Bio::Seq->new(-id=>'NOT3MUL', -seq=>'gatac');
 129 isa_ok($seq, 'Bio::PrimarySeqI');
 130 eval {
 131   Bio::Align::DNAStatistics->count_syn_sites($seq);
 132 };
 133 like($@, qr/not integral number of codons/);
 134
 135 # bug 2901
 136 $in = Bio::AlignIO->new(-file => test_input_file('bug2901.fa'),
 137                         -format => 'fasta');
 138
 139 $stats = Bio::Align::DNAStatistics->new(-verbose => 2);
 140 $aln = $in->next_aln();
 141 my $matrix;
 142 throws_ok {
 143 $matrix = $stats->distance(-align=>$aln,-method=>'Uncorrected');
 144 } qr/No distance calculated between seq3 and seq4/, "Warn if seqs don't overlap";
 145 $stats->verbose(-1);
 146 $matrix = $stats->distance(-align=>$aln,-method=>'Uncorrected');
 147 like($matrix->print_matrix, qr/-1/);
 148