use residue count for calcs instead of length, fixes bug #3185

bioperl · Apr 4, 2011 · 536fd50 · 536fd50
1 parent 5e4f151
commit 536fd50
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 13 deletions.
diff --git a/Bio/Tools/SeqStats.pm b/Bio/Tools/SeqStats.pm
@@ -571,16 +571,20 @@ sub get_mol_wt {
 	my $weight_lower_bound = 0;
 	my $weight_upper_bound = 0;
 	my $weight_table =  $Weights{$moltype};
-
+    my $total_res;
+
 	# compute weight of all the residues
 	foreach $element (keys %$rcount) {
 		$weight_lower_bound += $$rcount{$element} * $$weight_table{$element}->[0];
 		$weight_upper_bound += $$rcount{$element} * $$weight_table{$element}->[1];
+
+        # this tracks only the residues used for counting MW
+        $total_res += $$rcount{$element};
 	}
 	if ($moltype =~ /protein/) {
-    	# remove H2O during peptide bond formation.
-    	$weight_lower_bound -= $water * ($seqobj->length - 1);
-    	$weight_upper_bound -= $water * ($seqobj->length - 1);
+        # remove H2O during peptide bond formation.
+    	$weight_lower_bound -= $water * ($total_res - 1);
+    	$weight_upper_bound -= $water * ($total_res - 1);
 	} else {
     	# Correction because phosphate of 5' residue has additional OH and
     	# sugar ring of 3' residue has additional H

diff --git a/t/SeqTools/SeqStats.t b/t/SeqTools/SeqStats.t
@@ -7,7 +7,7 @@ BEGIN {
     use lib '.';
     use Bio::Root::Test;
 
-    test_begin(-tests => 44);
+    test_begin(-tests => 47);
 
 	use_ok('Bio::SeqIO');
 	use_ok('Bio::Tools::SeqStats');
@@ -49,7 +49,6 @@ 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);
@@ -62,7 +61,6 @@ $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
@@ -71,6 +69,8 @@ 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);
@@ -83,21 +83,33 @@ $wt = Bio::Tools::SeqStats->get_mol_wt($seqobj);
 is int $$wt[0], 2896;
 is int $$wt[1], 2896;
 
-$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;
+# Issue 3185: https://redmine.open-bio.org/issues/3185
 
-ok $seqobj = Bio::PrimarySeq->new(-seq=>'TGCCGTGTGTGCTGCTGCT', -alphabet=>'rna');
+$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 &round($$wt[0]), 6104 ;
+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
 #