Handle mixed alleles in get_single_allele_from_reads instead of crashing

Previously raised ValueError if reads had different alleles. Now uses the
most common allele and filters to only matching reads, with a warning.
Return signature changed to (allele, filtered_reads) tuple.

This handles edge cases from indel representation ambiguity or subclonal
mixtures without crashing. Callers using the old single-return-value
signature are updated (variant_sequence_creator.py).

Fixes #152

Author: iskandr

isovar/__init__.py

@@ -10,7 +10,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-__version__ = "1.4.13"
+__version__ = "1.4.14"
 
 
 from .allele_read import AlleleRead

isovar/allele_read_helpers.py

@@ -15,7 +15,7 @@
 AlleleRead objects.
 """
 
-from collections import defaultdict
+from collections import Counter, defaultdict
 
 from .common import groupby
 from .logging import get_logger
@@ -34,19 +34,38 @@ def group_reads_by_allele(allele_reads):
 
 def get_single_allele_from_reads(allele_reads):
     """
-    Given a sequence of AlleleRead objects, which are expected to all have
-    the same allele, return that allele.
+    Given a sequence of AlleleRead objects which are expected to all have
+    the same allele, return that allele and the reads supporting it.
+
+    If reads carry different alleles (e.g. from indel representation
+    ambiguity or subclonal mixtures), the most common allele is returned
+    along with only the reads that match it. A warning is logged about
+    the discarded minority alleles.
+
+    Returns (allele_string, list_of_matching_reads)
     """
     allele_reads = list(allele_reads)
 
     if len(allele_reads) == 0:
         raise ValueError("Expected non-empty list of AlleleRead objects")
 
-    seq = allele_reads[0].allele
-    if any(read.allele != seq for read in allele_reads):
-        raise ValueError("Expected all AlleleRead objects to have same allele '%s', got %s" % (
-            seq, allele_reads))
-    return seq
+    allele_counts = Counter(read.allele for read in allele_reads)
+    most_common_allele = allele_counts.most_common(1)[0][0]
+
+    if len(allele_counts) > 1:
+        logger.warning(
+            "Expected all reads to have same allele but found %d distinct "
+            "alleles: %s. Using most common allele '%s' (%d/%d reads).",
+            len(allele_counts),
+            dict(allele_counts),
+            most_common_allele,
+            allele_counts[most_common_allele],
+            len(allele_reads))
+        allele_reads = [
+            r for r in allele_reads if r.allele == most_common_allele
+        ]
+
+    return most_common_allele, allele_reads
 
 
 def group_unique_sequences(

isovar/variant_sequence_creator.py

@@ -89,7 +89,7 @@ def reads_to_variant_sequences(
         if len(variant_reads) == 0:
             return []
 
-        alt_seq = get_single_allele_from_reads(variant_reads)
+        alt_seq, variant_reads = get_single_allele_from_reads(variant_reads)
 
         # the number of context nucleotides on either side of the variant
         # is half the desired length (minus the number of variant nucleotides)

tests/test_read_helpers.py

@@ -11,10 +11,12 @@
 # limitations under the License.
 
 from isovar import ReadCollector
-from isovar.allele_read_helpers import  group_unique_sequences
+from isovar.allele_read import AlleleRead
+from isovar.allele_read_helpers import group_unique_sequences, get_single_allele_from_reads
 
 from varcode import Variant
 
+from .common import eq_
 from .testing_helpers import load_bam
 
 
@@ -45,3 +47,32 @@ def test_group_unique_sequences():
     # there are some redundant reads, so we expect that the number of
     # unique entries should be less than the total read partitions
     assert len(variant_reads) > len(groups)
+
+
+def test_get_single_allele_from_reads_uniform():
+    reads = [
+        AlleleRead(prefix="AA", allele="G", suffix="TT", name="r1"),
+        AlleleRead(prefix="AA", allele="G", suffix="TT", name="r2"),
+    ]
+    allele, filtered_reads = get_single_allele_from_reads(reads)
+    eq_(allele, "G")
+    eq_(len(filtered_reads), 2)
+
+
+def test_get_single_allele_from_reads_mixed_uses_most_common():
+    reads = [
+        AlleleRead(prefix="AA", allele="G", suffix="TT", name="r1"),
+        AlleleRead(prefix="AA", allele="G", suffix="TT", name="r2"),
+        AlleleRead(prefix="AA", allele="G", suffix="TT", name="r3"),
+        AlleleRead(prefix="AA", allele="C", suffix="TT", name="r4"),
+    ]
+    allele, filtered_reads = get_single_allele_from_reads(reads)
+    eq_(allele, "G")
+    eq_(len(filtered_reads), 3)
+    assert all(r.allele == "G" for r in filtered_reads)
+
+
+def test_get_single_allele_from_reads_empty_raises():
+    import pytest
+    with pytest.raises(ValueError):
+        get_single_allele_from_reads([])