t/Align/AlignStats.t

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