# -*-Perl-*- Test Harness script for Bioperl
# $Id: SeqStats.t 11525 2007-06-27 10:16:38Z sendu $
use strict;
BEGIN {
use lib '.';
use Bio::Root::Test;
test_begin(-tests => 47);
use_ok('Bio::SeqIO');
use_ok('Bio::Tools::SeqStats');
}
my ($seqobj, $count, $seqobj_stats, $wt);
my $str = Bio::SeqIO->new(-file=> test_input_file('multifa.seq'), '-format' => 'Fasta');
$seqobj = $str->next_seq();
ok defined $seqobj, 'new Bio::Root::IO object';
my $seq_stats = Bio::Tools::SeqStats->new('-seq' => $seqobj);
ok defined $seq_stats && $seq_stats, 'new Bio::Tools:SeqStats object';
# eg for DNA sequence
my $hash_ref = $seq_stats->count_monomers();
is ( $hash_ref->{'A'}, 80 , 'count_monomers()');
$hash_ref = $seq_stats-> count_codons();
ok defined $hash_ref && $hash_ref , 'count_codons()';
my $weight = $seq_stats->get_mol_wt();
ok defined $weight && $weight , 'get_mol_wt()' ;
$seqobj = Bio::PrimarySeq->new(-seq=>'ACTGTGGCGTCAACTG',
-alphabet=>'dna', -id=>'test');
ok $seqobj_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj);
isa_ok $seqobj_stats, 'Bio::Tools::SeqStats';
$count = $seqobj_stats->count_monomers(); # for DNA sequence
is $count->{'A'}, 3;
is $count->{'C'}, 4;
is $count->{'G'}, 5;
is $count->{'T'}, 4;
$count = $seqobj_stats->count_codons();
is $count->{'ACT'}, 2;
is $count->{'GTG'}, 1;
is $count->{'GCG'}, 1;
is $count->{'TCA'}, 1;
$seqobj = Bio::PrimarySeq->new(-seq=>'ACTACTTCA', -alphabet=>'dna',
-id=>'test');
$seqobj_stats = Bio::Tools::SeqStats->new('-seq' => $seqobj);
$wt = $seqobj_stats->get_mol_wt(); # for DNA sequence
is &round($$wt[0]), 2738 ;
$seqobj = Bio::PrimarySeq->new(-seq=>'ACXACNNCA',
-alphabet=>'dna', -id=>'test');
$wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
is &round($$wt[0]), 2693;
is &round($$wt[1]), 2813;
$seqobj = Bio::PrimarySeq->new(-seq=>'ACTGTGGCGTCAACTG',
-alphabet=>'dna', -id=>'test');
$count = Bio::Tools::SeqStats->count_monomers($seqobj); # for DNA sequence
is $count->{'A'}, 3;
is $count->{'C'}, 4;
is $count->{'G'}, 5;
is $count->{'T'}, 4;
# protein
$seqobj = Bio::PrimarySeq->new(-seq=>'MQSERGITIDISLWKFETSKYYVT',
-alphabet=>'protein', -id=>'test');
$seqobj_stats = Bio::Tools::SeqStats->new('-seq' => $seqobj);
$count = $seqobj_stats->count_monomers(); # for amino sequence
is $$count{'M'}, 1;
is $$count{'I'}, 3;
is $$count{'Y'}, 2;
is $$count{'T'}, 3;
$wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
is int $$wt[0], 2896;
is int $$wt[1], 2896;
# Issue 3185: https://redmine.open-bio.org/issues/3185
$seqobj = Bio::PrimarySeq->new(-seq=>'MQSERGITIDISLWKFETSKYYVT*',
-alphabet=>'protein', -id=>'test');
is($seqobj->seq, 'MQSERGITIDISLWKFETSKYYVT*');
$seqobj_stats = Bio::Tools::SeqStats->new('-seq' => $seqobj);
$wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
is int $$wt[0], 2896;
is int $$wt[1], 2896;
# selenocysteine
ok $seqobj = Bio::PrimarySeq->new(-seq=>'MQSERGITIDISLWKFETSKYYVT',
-alphabet=>'protein');
$wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
is &round($$wt[0]), 2896 ;
# RNA
$seqobj = Bio::PrimarySeq->new(-seq=>'UYXUYNNYU', -alphabet=>'rna');
$wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
is &round($$wt[0]), 2768;
is &round($$wt[1]), 2891;
ok $seqobj = Bio::PrimarySeq->new(-seq=>'TGCCGTGTGTGCTGCTGCT', -alphabet=>'rna');
$wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
is &round($$wt[0]), 6104 ;
#
# hydropathicity aka "gravy" score
#
# normal seq (should succeed)
ok $seqobj = Bio::PrimarySeq->new(-seq=>'MSFVLVAPDMLATAAADVVQIGSAVSAGS',
-alphabet=>'protein');
my $gravy = Bio::Tools::SeqStats->hydropathicity($seqobj);
is int($gravy*1000), 1224; # check to nearest 0.1%
# ambiguous sequence (should fail)
ok $seqobj = Bio::PrimarySeq->new(-seq=>'XXXB**BS', -alphabet=>'protein');
eval { Bio::Tools::SeqStats->hydropathicity($seqobj) };
like $@, qr/ambiguous amino acids/i;
# empty sequence (should fail)
ok $seqobj = Bio::PrimarySeq->new(-seq=>'', -alphabet=>'protein');
eval { Bio::Tools::SeqStats->hydropathicity($seqobj) };
like $@, qr/hydropathicity not defined/i;
# DNA sequence (should fail)
ok $seqobj = Bio::PrimarySeq->new(-seq=>'GATTACA', -alphabet=>'dna');
eval { Bio::Tools::SeqStats->hydropathicity($seqobj) };
like $@, qr/only meaningful for protein/;
#
# Extra functions
#
# perl does not have an explicit rounding function
sub round { return int ((shift @_) + 0.5 ) }