Repository: acorg/dark-matter

Issue:
SAM referenceInsertion should record what reads cause reference insertions
The code currently records all reference insertions but doesn't keep track of which reads caused insertions. That makes it impossible to ignore insertions according to which read created it.

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d7f91b4..add466a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,9 @@
+## 3.0.23 June 22, 2018
+
+Changed the way reference sequence insertions are stored in a
+`dark.sam.PaddedSAM` instance to make it possible to tell which query
+sequences caused reference insertions.
+
 ## 3.0.22 June 21, 2018
 
 Made `dark/sam.py` properly deal with secondary alignments that are missing
diff --git a/dark/__init__.py b/dark/__init__.py
index db6691a..e3d166b 100644
--- a/dark/__init__.py
+++ b/dark/__init__.py
@@ -7,4 +7,4 @@ if sys.version_info < (2, 7):
 # will not be found by the version() function in ../setup.py
 #
 # Remember to update ../CHANGELOG.md describing what's new in each version.
-__version__ = '3.0.22'
+__version__ = '3.0.23'
diff --git a/dark/sam.py b/dark/sam.py
index f0ffa19..2165ea6 100644
--- a/dark/sam.py
+++ b/dark/sam.py
@@ -29,12 +29,12 @@ class PaddedSAM(object):
     """
     def __init__(self, filename):
         self.samfile = AlignmentFile(filename)
-        # self.referenceInsertions will be keyed by offset into the reference
-        # sequence. The inserted bases would need to begin at this offset. The
-        # value will be a Counter whose keys are the nucleotides proposed for
-        # insertion, with a value indicating how many times the nucleotide was
-        # proposed for insertion at that offset.
-        self.referenceInsertions = defaultdict(Counter)
+        # self.referenceInsertions will be keyed by query id (the query
+        # that would cause a reference insertion). The values will be lists
+        # of 2-tuples, with each 2-tuple containing an offset into the
+        # reference sequence and the C{str} of nucleotide that would be
+        # inserted starting at that offset.
+        self.referenceInsertions = defaultdict(list)
 
     def close(self):
         """
@@ -182,6 +182,16 @@ class PaddedSAM(object):
                 if rcSuffix:
                     read.query_name += rcSuffix
 
+            # Adjust the query id if it's a duplicate and we're not
+            # allowing duplicates.
+            if allowDuplicateIds:
+                queryId = read.query_name
+            else:
+                count = idCount[read.query_name]
+                idCount[read.query_name] += 1
+                queryId = read.query_name + (
+                    '' if count == 0 else '/%d' % count)
+
             referenceStart = read.reference_start
             atStart = True
             queryIndex = 0
@@ -205,9 +215,9 @@ class PaddedSAM(object):
                     # query but record what would have been inserted into the
                     # reference.
                     atStart = False
-                    for i in range(length):
-                        self.referenceInsertions[referenceIndex + i][
-                            query[queryIndex + i]] += 1
+                    self.referenceInsertions[queryId].append(
+                        (referenceIndex,
+                         query[queryIndex:queryIndex + length]))
                 elif operation == CDEL:
                     # Delete from the reference. Some bases from the reference
                     # would need to be deleted to continue the match. So we put
@@ -293,14 +303,7 @@ class PaddedSAM(object):
                 padChar * (referenceLength -
                            (referenceStart + len(alignedSequence))))
 
-            if allowDuplicateIds:
-                suffix = ''
-            else:
-                count = idCount[read.query_name]
-                idCount[read.query_name] += 1
-                suffix = '' if count == 0 else '/%d' % count
-
-            yield Read('%s%s' % (read.query_name, suffix), paddedSequence)
+            yield Read(queryId, paddedSequence)
 
 
 @contextmanager


Test patch:
diff --git a/test/test_sam.py b/test/test_sam.py
index e31f801..f8fd7e9 100644
--- a/test/test_sam.py
+++ b/test/test_sam.py
@@ -347,8 +347,32 @@ class TestPaddedSAM(TestCase):
             self.assertEqual(Read('query1', '-TCGG-----'), read)
             self.assertEqual(
                 {
-                    3: {'T': 1},
-                    4: {'A': 1},
+                    'query1': [(3, 'TA')],
+                },
+                ps.referenceInsertions)
+            ps.close()
+
+    def testPrimaryAndSecondaryReferenceInsertion(self):
+        """
+        A primary and secondary insertion into the reference (of the same
+        query) must result in the expected padded sequences and the expected
+        value in the referenceInsertions dictionary.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            'query1 0 ref1 2 60 2M2I2M * 0 0 TCTAGG ZZZZZZ',
+            'query1 256 ref1 4 60 2M3I1M * 0 0 * *',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            ps = PaddedSAM(filename)
+            (read1, read2) = list(ps.queries())
+            self.assertEqual(Read('query1', '-TCGG-----'), read1)
+            self.assertEqual(Read('query1/1', '---TCG----'), read2)
+            self.assertEqual(
+                {
+                    'query1': [(3, 'TA')],
+                    'query1/1': [(5, 'TAG')],
                 },
                 ps.referenceInsertions)
             ps.close()