Fix boundary read evidence counts

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.21"
+__version__ = "1.4.22"
 
 
 from .allele_read import AlleleRead

isovar/read_collector.py

@@ -152,6 +152,9 @@ def locus_read_from_pysam_aligned_segment(
             if ref_pos is not None
         }
 
+        aligned_subsequence_start = pysam_aligned_segment.query_alignment_start
+        aligned_subsequence_end = pysam_aligned_segment.query_alignment_end
+
         reference_interval_size = base0_end_exclusive - base0_start_inclusive
         if reference_interval_size < 0:
             raise ValueError("Unexpected interval start after interval end")
@@ -194,16 +197,22 @@ def locus_read_from_pysam_aligned_segment(
             read_base0_before_insertion = reference_position_to_read_position.get(
                 reference_position_before_insertion
             )
-            if read_base0_before_insertion is None:
-                return None
-
             read_base0_after_insertion = reference_position_to_read_position.get(
                 reference_position_after_insertion
             )
-            if read_base0_after_insertion is None:
-                return None
 
-            if read_base0_after_insertion - read_base0_before_insertion == 1:
+            if (
+                read_base0_before_insertion is None
+                and read_base0_after_insertion is None
+            ):
+                return None
+            elif read_base0_before_insertion is None:
+                read_base0_start_inclusive = aligned_subsequence_start
+                read_base0_end_exclusive = read_base0_after_insertion
+            elif read_base0_after_insertion is None:
+                read_base0_start_inclusive = read_base0_before_insertion + 1
+                read_base0_end_exclusive = aligned_subsequence_end
+            elif read_base0_after_insertion - read_base0_before_insertion == 1:
                 read_base0_start_inclusive = read_base0_end_exclusive = (
                     read_base0_before_insertion + 1
                 )
@@ -269,8 +278,6 @@ def locus_read_from_pysam_aligned_segment(
             # than the sequence, qualities, and alignment positions
             # we've extracted, so slice through those to get rid of
             # soft-clipped ends of the read
-            aligned_subsequence_start = pysam_aligned_segment.query_alignment_start
-            aligned_subsequence_end = pysam_aligned_segment.query_alignment_end
             sequence = sequence[aligned_subsequence_start:aligned_subsequence_end]
             base0_reference_positions = base0_reference_positions[
                 aligned_subsequence_start:aligned_subsequence_end

tests/test_cli.py

@@ -77,7 +77,7 @@ def test_cli_allele_reads():
     df = run_cli_fn(isovar_allele_reads, return_dataframe=True)
     assert set(["prefix", "allele", "suffix", "name", "sequence", "gene"]).issubset(df.columns)
     assert "ref_reads" not in df.columns
-    assert len(df) == 293
+    assert len(df) == 294
 
 
 def test_cli_reference_contexts():

tests/test_locus_reads.py

@@ -235,13 +235,9 @@ def test_locus_reads_dataframe():
         for line in f:
             if line.startswith("HWI"):
                 n_reads_expected += 1
-    # we know from inspecting the file that *one* of the reads overlapping this
-    # variant has a CIGAR string of N at the location before and thus we'll
-    # be missing that read.
-    #
-    # TODO: figure out what to do when the variant nucleotide is at the start or
-    # end of an exon, since that won't have mapping positions on both its left
-    # and right
+    # We know from inspecting the file that one fetched read is spliced across
+    # this locus and therefore does not actually overlap the queried interval.
+    # The widened fetch window pulls it in, and get_overlap() correctly drops it.
     n_reads_expected -= 1
 
     print("Found %d sequences in %s" % (n_reads_expected, sam_path_single_variant))
@@ -410,3 +406,51 @@ def test_locus_read_insertion_locus_with_insertion_in_read():
     assert result is not None, "Read with insertion should be returned"
     eq_(result.read_base0_start_inclusive, 4)
     eq_(result.read_base0_end_exclusive, 5)
+
+
+def test_locus_read_insertion_locus_with_insertion_at_start_of_read():
+    """
+    An insertion-supporting read can begin at the insertion locus and still
+    provide a non-empty inserted interval.
+
+    Regression test for GitHub issue #49.
+    """
+    pysam_read = make_pysam_read(
+        seq="GT",
+        cigar="1I1M",
+        mdtag="1",
+        reference_start=3,
+    )
+    read_collector = ReadCollector()
+    result = read_collector.locus_read_from_pysam_aligned_segment(
+        pysam_read,
+        base0_start_inclusive=3,
+        base0_end_exclusive=3,
+    )
+    assert result is not None, "Insertion at the start of a read should be kept"
+    eq_(result.read_base0_start_inclusive, 0)
+    eq_(result.read_base0_end_exclusive, 1)
+
+
+def test_locus_read_insertion_locus_with_insertion_at_end_of_read():
+    """
+    An insertion-supporting read can end at the insertion locus and still
+    provide a non-empty inserted interval.
+
+    Regression test for GitHub issue #49.
+    """
+    pysam_read = make_pysam_read(
+        seq="AG",
+        cigar="1M1I",
+        mdtag="1",
+        reference_start=2,
+    )
+    read_collector = ReadCollector()
+    result = read_collector.locus_read_from_pysam_aligned_segment(
+        pysam_read,
+        base0_start_inclusive=3,
+        base0_end_exclusive=3,
+    )
+    assert result is not None, "Insertion at the end of a read should be kept"
+    eq_(result.read_base0_start_inclusive, 1)
+    eq_(result.read_base0_end_exclusive, 2)

tests/test_variant_reads_with_dummy_samfile.py

@@ -133,3 +133,164 @@ def test_partitioned_read_sequences_deletion():
         allele="",
         suffix="G")
     eq_(variant_read, expected)
+
+
+def test_read_evidence_counts_snv_reads_at_read_boundaries():
+    """
+    Reads whose SNV locus is the first or last aligned base should still be
+    counted.
+
+    Regression test for GitHub issue #55.
+    """
+    chromosome = "1"
+    variant = Variant(
+        chromosome,
+        4,
+        "T",
+        "G",
+        grch38,
+        normalize_contig_names=False)
+
+    reads = [
+        make_pysam_read(
+            seq="GAAA",
+            cigar="4M",
+            mdtag="0T3",
+            name="alt-start",
+            reference_start=3),
+        make_pysam_read(
+            seq="TAAA",
+            cigar="4M",
+            mdtag="4",
+            name="ref-start",
+            reference_start=3),
+        make_pysam_read(
+            seq="AAAG",
+            cigar="4M",
+            mdtag="3T0",
+            name="alt-end",
+            reference_start=0),
+        make_pysam_read(
+            seq="AAAT",
+            cigar="4M",
+            mdtag="4",
+            name="ref-end",
+            reference_start=0),
+    ]
+
+    read_creator = ReadCollector()
+    read_evidence = read_creator.read_evidence_for_variant(
+        variant=variant,
+        alignment_file=MockAlignmentFile(references=(chromosome,), reads=reads))
+
+    eq_(read_evidence.alt_read_names, {"alt-start", "alt-end"})
+    eq_(read_evidence.ref_read_names, {"ref-start", "ref-end"})
+
+
+def test_read_evidence_counts_insertion_reads_at_read_boundaries():
+    """
+    Insertion-supporting reads that begin or end at the insertion locus should
+    be counted as alt reads instead of being dropped.
+
+    Regression test for GitHub issue #49. Also serves as an indel-count
+    regression for GitHub issue #23.
+    """
+    chromosome = "1"
+    variant = Variant(
+        chromosome,
+        3,
+        "A",
+        "AG",
+        grch38,
+        normalize_contig_names=False)
+
+    reads = [
+        make_pysam_read(
+            seq="GT",
+            cigar="1I1M",
+            mdtag="1",
+            name="alt-start",
+            reference_start=3),
+        make_pysam_read(
+            seq="T",
+            cigar="1M",
+            mdtag="1",
+            name="ref-start",
+            reference_start=3),
+        make_pysam_read(
+            seq="AG",
+            cigar="1M1I",
+            mdtag="1",
+            name="alt-end",
+            reference_start=2),
+        make_pysam_read(
+            seq="A",
+            cigar="1M",
+            mdtag="1",
+            name="ref-end",
+            reference_start=2),
+    ]
+
+    read_creator = ReadCollector()
+    read_evidence = read_creator.read_evidence_for_variant(
+        variant=variant,
+        alignment_file=MockAlignmentFile(references=(chromosome,), reads=reads))
+
+    eq_(read_evidence.alt_read_names, {"alt-start", "alt-end"})
+    eq_(read_evidence.ref_read_names, {"ref-start", "ref-end"})
+
+
+def test_partitioned_read_sequences_snv_at_last_exonic_base_before_splice():
+    """
+    A variant on the last nucleotide of an exon should still be recovered from
+    a spliced read.
+
+    Regression test for GitHub issue #24.
+    """
+    chromosome = "1"
+    variant = Variant(
+        chromosome,
+        4,
+        "T",
+        "G",
+        grch38,
+        normalize_contig_names=False)
+    read = make_pysam_read(
+        seq="ACCGTGGA",
+        cigar="4M100N4M",
+        name="splice-before",
+        reference_start=0)
+    read_creator = ReadCollector()
+    variant_reads = read_creator.allele_reads_supporting_variant(
+        variant=variant,
+        alignment_file=MockAlignmentFile(references=(chromosome,), reads=[read]))
+    assert len(variant_reads) == 1
+    eq_(variant_reads[0].allele, "G")
+
+
+def test_partitioned_read_sequences_snv_at_first_exonic_base_after_splice():
+    """
+    A variant on the first nucleotide of an exon should still be recovered from
+    a spliced read.
+
+    Regression test for GitHub issue #24.
+    """
+    chromosome = "1"
+    variant = Variant(
+        chromosome,
+        105,
+        "T",
+        "G",
+        grch38,
+        normalize_contig_names=False)
+    read = make_pysam_read(
+        seq="ACCTGGGA",
+        cigar="4M100N4M",
+        name="splice-after",
+        reference_start=0)
+    read_creator = ReadCollector()
+    variant_reads = read_creator.allele_reads_supporting_variant(
+        variant=variant,
+        alignment_file=MockAlignmentFile(references=(chromosome,), reads=[read]))
+    assert len(variant_reads) == 1
+    eq_(variant_reads[0].allele, "G")