Repository: acorg/dark-matter

Issue:
Change SAM filtering to allow for CIGAR strings that indicate hard clipping but that have sequences that are not clipped
It can happen that an aligner (`bwa mem -a` does this) can generate SAM lines like these two:

```
K00234:90:HTWVHBBXX:6:1101:1844:10493 0 Chimp-D00220	364 0 3S78M * 0 0
TTTTGGTTATCGCTGGATGTGTCTGCGGCGTTTTATCATCTTCCTCTTCATCCTGCTGCTATGCCTCATCTTATTGTTGGT
AAFFFJAJJJJJFJJJ7JJJJFJFJJFJFJAFJJJJAAF<FFFAF<JJ-FF<JF-AF7F<AJF--<-F7-AA-7JFJJFJ<
NM:i:1 MD:Z:69C8 AS:i:73 XS:i:73

K00234:90:HTWVHBBXX:6:1101:1844:10493 256 D-AM494716 364 0 3H78M * 0 0
* *
NM:i:1 MD:Z:69C8 AS:i:73
```

those are actually just 2 lines in the SAM file, giving alignments for the same read. The issue here is that the first alignment has a CIGAR string of `3S78M` (3 soft clipped nts, 78 matches) and because it's a soft clip the bases are still in the query string. On the next line we have `*` and `*` for the query and quality (because this is a secondary match (flag = 256)) and a CIGAR string of `3H78M` (3 hard clipped, 78 match). Normally with hard clipping the bases are removed from the sequence, but in this case because the query is first given with a soft clip the bases are retained. So code that assumes hard -clipped bases will be removed will break.

You might think `bwa mem -a -H` would cause hard clipping to always be done (perhaps by repeating the (clipped) query sequence in secondary matches) but that's not what happens - bwa reverts to the non- `-a` option behaviour.

So our code in `dark/sam.py` will need to be smarter. It should figure out in advance if the query has actually been hard clipped and use that to decide whether a hard clip indicator in the CIGAR string should consume query bases or not.

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d82fc9b..0e884c1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,8 @@
+## 3.0.39 October 10, 2018
+
+The fix to solve [#630](https://github.com/acorg/dark-matter/issues/630)
+was insufficient. That's fixed in this release, hopefully!
+
 ## 3.0.38 October 5, 2018
 
 Fixed [#630](https://github.com/acorg/dark-matter/issues/630) to deal with
diff --git a/dark/__init__.py b/dark/__init__.py
index 7513e6e..20729ae 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.38'
+__version__ = '3.0.39'
diff --git a/dark/sam.py b/dark/sam.py
index 0167fdc..a4f2c1a 100644
--- a/dark/sam.py
+++ b/dark/sam.py
@@ -65,40 +65,56 @@ def samReferencesToStr(filenameOrSamfile, indent=0):
         return _references(sam)
 
 
-def _hardClip(sequence, cigartuples):
+def _hardClip(sequence, quality, cigartuples):
     """
     Hard clip (if necessary) a sequence.
 
     @param sequence: A C{str} nucleotide sequence.
-    @param cigartuples: An iterable of (operation, length) tuples, specifying
-        matching as per the SAM specification.
-    @return: A hard-clipped C{str} sequence if hard-clipping is indicated by
-        the CIGAR operations and has not already been performed (as indicated
-        by the lengths of the sequence and the sum of the CIGAR operation
-        lengths).
+    @param quality: A C{str} quality string or a C{list} of C{int} quality
+        values as returned by pysam.
+    @param cigartuples: An iterable of (operation, length) tuples, detailing
+        the alignment, as per the SAM specification.
+    @return: A 3-tuple consisting of
+            1) a hard-clipped C{str} sequence if hard-clipping is indicated by
+               the CIGAR operations.
+            2) a hard-clipped quality C{str} or C{list} (depending on what
+               type we were passed) if hard-clipping is indicated by the CIGAR
+               operations.
+            3) a Boolean, C{True} if hard clipping was performed by this
+               function or C{False} if the hard clipping had already been
+               done.
     """
     hardClipCount = cigarLength = 0
     for (operation, length) in cigartuples:
         hardClipCount += operation == CHARD_CLIP
-        cigarLength += length
+        cigarLength += length if operation in _CONSUMES_QUERY else 0
 
     sequenceLength = len(sequence)
+    assert sequenceLength == len(quality)
     clipLeft = clipRight = 0
+    clippedSequence = sequence
+    clippedQuality = quality
+
+    if sequenceLength > cigarLength:
+        alreadyClipped = False
+    else:
+        assert sequenceLength == cigarLength
+        alreadyClipped = True
 
     if hardClipCount == 0:
         pass
     elif hardClipCount == 1:
         # Hard clip either at the start or the end.
         if cigartuples[0][0] == CHARD_CLIP:
-            clipLeft = cigartuples[0][1]
-            if sequenceLength == cigarLength:
-                # The LHS hard clipping has not been done.
-                sequence = sequence[clipLeft:]
+            if not alreadyClipped:
+                clipLeft = cigartuples[0][1]
+                clippedSequence = sequence[clipLeft:]
+                clippedQuality = quality[clipLeft:]
         elif cigartuples[-1][0] == CHARD_CLIP:
-            clipRight = cigartuples[-1][1]
-            if sequenceLength == cigarLength:
-                # The RHS hard clipping has not been done.
-                sequence = sequence[:-clipRight]
+            if not alreadyClipped:
+                clipRight = cigartuples[-1][1]
+                clippedSequence = sequence[:-clipRight]
+                clippedQuality = quality[:-clipRight]
         else:
             raise ValueError(
                 'Invalid CIGAR tuples (%s) contains hard-clipping operation '
@@ -107,19 +123,33 @@ def _hardClip(sequence, cigartuples):
     elif hardClipCount == 2:
         # Hard clip at both the start and end.
         assert cigartuples[0][0] == cigartuples[-1][0] == CHARD_CLIP
-        clipLeft, clipRight = cigartuples[0][1], cigartuples[-1][1]
-        if sequenceLength == cigarLength:
-            # The hard clipping has not been done.
-            sequence = sequence[clipLeft:-clipRight]
+        if not alreadyClipped:
+            clipLeft, clipRight = cigartuples[0][1], cigartuples[-1][1]
+            clippedSequence = sequence[clipLeft:-clipRight]
+            clippedQuality = quality[clipLeft:-clipRight]
     else:
         raise ValueError(
             'Invalid CIGAR tuples (%s) specifies hard-clipping %d times (2 '
             'is the maximum).' % (cigartuples, hardClipCount))
 
-    assert len(sequence) + clipLeft + clipRight == cigarLength, (
-        '%d + %d + %d != %d' % (len(sequence), clipLeft, clipRight,
-                                cigarLength))
-    return sequence
+    weClipped = bool(clipLeft or clipRight)
+
+    if weClipped:
+        assert not alreadyClipped
+        if len(clippedSequence) + clipLeft + clipRight != sequenceLength:
+            raise ValueError(
+                'Sequence %r (length %d) clipped to %r (length %d), but the '
+                'difference between these two lengths (%d) is not equal to '
+                'the sum (%d) of the left and right clip lengths (%d and %d '
+                'respectively). CIGAR tuples: %s' %
+                (sequence, len(sequence),
+                 clippedSequence, len(clippedSequence),
+                 abs(len(sequence) - len(clippedSequence)),
+                 clipLeft + clipRight, clipLeft, clipRight, cigartuples))
+    else:
+        assert len(clippedSequence) == len(clippedQuality) == sequenceLength
+
+    return clippedSequence, clippedQuality, weClipped
 
 
 class SAMFilter(object):
@@ -280,6 +310,7 @@ class SAMFilter(object):
         if storeQueryIds:
             self.queryIds = queryIds = set()
 
+        lastAlignment = None
         count = 0
         with samfile(self.filename) as samAlignment:
             for count, alignment in enumerate(samAlignment.fetch(), start=1):
@@ -296,6 +327,32 @@ class SAMFilter(object):
                                 (maxScore is not None and score > maxScore)):
                             continue
 
+                # Secondary and supplementary alignments may have a '*'
+                # (pysam returns this as None) SEQ field, indicating that
+                # the previous sequence should be used. This is best
+                # practice according to section 2.5.2 of
+                # https://samtools.github.io/hts-specs/SAMv1.pdf So we use
+                # the last alignment if we get None as a query sequence.
+                if alignment.query_sequence is None:
+                    if lastAlignment is None:
+                        raise InvalidSAM(
+                            'pysam produced an alignment (number %d) with no '
+                            'query sequence without previously giving an '
+                            'alignment with a sequence.' % count)
+                    # Use the previous query sequence and quality.
+                    (alignment.query_sequence,
+                     alignment.query_qualities, _) = _hardClip(
+                         lastAlignment.query_sequence,
+                         lastAlignment.query_qualities,
+                         alignment.cigartuples)
+                else:
+                    lastAlignment = alignment
+                    (alignment.query_sequence,
+                     alignment.query_qualities, _) = _hardClip(
+                         alignment.query_sequence,
+                         alignment.query_qualities,
+                         alignment.cigartuples)
+
                 if ((filterRead is None or
                      filterRead(Read(alignment.query_name,
                                      alignment.query_sequence,
@@ -372,8 +429,8 @@ class PaddedSAM(object):
         self.referenceInsertions = defaultdict(list)
 
     def queries(self, rcSuffix='', rcNeeded=False, padChar='-',
-                queryInsertionChar='N', allowDuplicateIds=False,
-                addAlignment=False):
+                queryInsertionChar='N', unknownQualityChar='!',
+                allowDuplicateIds=False, addAlignment=False):
         """
         Produce padded (with gaps) queries according to the CIGAR string and
         reference sequence length for each matching query sequence.
@@ -395,6 +452,9 @@ class PaddedSAM(object):
             is inserted as a 'missing' query character (i.e., a base that can
             be assumed to have been lost due to an error) whose existence is
             necessary for the match to continue.
+        @param unknownQualityChar: The character to put into the quality
+            string when unknown bases are inserted in the query or the query
+            is padded on the left/right with gaps.
         @param allowDuplicateIds: If C{True}, repeated query ids (due to
             secondary or supplemental matches) will not have /1, /2, etc.
             appended to their ids. So repeated ids may appear in the yielded
@@ -417,38 +477,17 @@ class PaddedSAM(object):
         idCount = Counter()
 
         MATCH_OPERATIONS = {CMATCH, CEQUAL, CDIFF}
-        lastQuery = None
 
         for lineNumber, alignment in enumerate(
                 self.samFilter.alignments(), start=1):
 
             query = alignment.query_sequence
-
-            # Secondary (and presumably supplementary) alignments may have
-            # a '*' (None in pysam) SEQ field, indicating that the previous
-            # sequence should be used. This is best practice according to
-            # section 2.5.2 of https://samtools.github.io/hts-specs/SAMv1.pdf
-            if query is None:
-                if alignment.is_secondary or alignment.is_supplementary:
-                    if lastQuery is None:
-                        raise InvalidSAM(
-                            'Query line %d has an empty SEQ field, but no '
-                            'previous alignment is present.' % lineNumber)
-                    else:
-                        query = _hardClip(lastQuery, alignment.cigartuples)
-                else:
-                    raise InvalidSAM(
-                        'Query line %d has an empty SEQ field, but the '
-                        'alignment is not marked as secondary or '
-                        'supplementary.' % lineNumber)
-            else:
-                # Remember the last query here (before we potentially modify
-                # it due to it being reverse complimented for the alignment).
-                lastQuery = query
+            quality = ''.join(chr(q + 33) for q in alignment.query_qualities)
 
             if alignment.is_reverse:
                 if rcNeeded:
                     query = DNARead('id', query).reverseComplement().sequence
+                    quality = quality[::-1]
                 if rcSuffix:
                     alignment.query_name += rcSuffix
 
@@ -467,6 +506,7 @@ class PaddedSAM(object):
             queryIndex = 0
             referenceIndex = referenceStart
             alignedSequence = ''
+            alignedQuality = ''
 
             for operation, length in alignment.cigartuples:
 
@@ -477,6 +517,7 @@ class PaddedSAM(object):
                 if operation in MATCH_OPERATIONS:
                     atStart = False
                     alignedSequence += query[queryIndex:queryIndex + length]
+                    alignedQuality += quality[queryIndex:queryIndex + length]
                 elif operation == CINS:
                     # Insertion to the reference. This consumes query bases but
                     # we don't output them because the reference cannot be
@@ -494,6 +535,7 @@ class PaddedSAM(object):
                     # an insertion into the query to compensate.
                     atStart = False
                     alignedSequence += queryInsertionChar * length
+                    alignedQuality += unknownQualityChar * length
                 elif operation == CREF_SKIP:
                     # Skipped reference. Opens a gap in the query. For
                     # mRNA-to-genome alignment, an N operation represents an
@@ -502,6 +544,7 @@ class PaddedSAM(object):
                     # to occur.
                     atStart = False
                     alignedSequence += queryInsertionChar * length
+                    alignedQuality += unknownQualityChar * length
                 elif operation == CSOFT_CLIP:
                     # Bases in the query that are not part of the match. We
                     # remove these from the query if they protrude before the
@@ -514,13 +557,17 @@ class PaddedSAM(object):
                         unwantedLeft = length - referenceStart
                         if unwantedLeft > 0:
                             # The query protrudes left. Copy its right part.
-                            alignedSequence += query[queryIndex + unwantedLeft:
-                                                     queryIndex + length]
+                            alignedSequence += query[
+                                queryIndex + unwantedLeft:queryIndex + length]
+                            alignedQuality += quality[
+                                queryIndex + unwantedLeft:queryIndex + length]
                             referenceStart = 0
                         else:
                             referenceStart -= length
                             alignedSequence += query[
                                 queryIndex:queryIndex + length]
+                            alignedQuality += quality[
+                                queryIndex:queryIndex + length]
                     else:
                         unwantedRight = (
                             (referenceStart + len(alignedSequence) + length) -
@@ -530,9 +577,13 @@ class PaddedSAM(object):
                             # The query protrudes right. Copy its left part.
                             alignedSequence += query[
                                 queryIndex:queryIndex + length - unwantedRight]
+                            alignedQuality += quality[
+                                queryIndex:queryIndex + length - unwantedRight]
                         else:
                             alignedSequence += query[
                                 queryIndex:queryIndex + length]
+                            alignedQuality += quality[
+                                queryIndex:queryIndex + length]
                 elif operation == CHARD_CLIP:
                     # Some bases have been completely removed from the query.
                     # This (H) can only be present as the first and/or last
@@ -555,8 +606,8 @@ class PaddedSAM(object):
             if queryIndex != len(query):
                 # Oops, we did not consume the entire query.
                 raise ValueError(
-                    'Query %s not fully consumed when parsing CIGAR string. '
-                    'Query %s (len %d), final query index %d, CIGAR: %r' %
+                    'Query %r not fully consumed when parsing CIGAR string. '
+                    'Query %r (len %d), final query index %d, CIGAR: %r' %
                     (alignment.query_name, query, len(query), queryIndex,
                      alignment.cigartuples))
 
@@ -567,13 +618,17 @@ class PaddedSAM(object):
 
             # Put gap characters before and after the aligned sequence so that
             # it is offset properly and matches the length of the reference.
-            paddedSequence = (
-                (padChar * referenceStart) +
-                alignedSequence +
-                padChar * (referenceLength -
-                           (referenceStart + len(alignedSequence))))
+            padRightLength = (referenceLength -
+                              (referenceStart + len(alignedSequence)))
+            paddedSequence = (padChar * referenceStart +
+                              alignedSequence +
+                              padChar * padRightLength)
+            paddedQuality = (unknownQualityChar * referenceStart +
+                             alignedQuality +
+                             unknownQualityChar * padRightLength)
+
+            read = Read(queryId, paddedSequence, paddedQuality)
 
-            read = Read(queryId, paddedSequence)
             if addAlignment:
                 read.alignment = alignment
 


Test patch:
diff --git a/test/test_sam.py b/test/test_sam.py
index 35acba8..4abf556 100644
--- a/test/test_sam.py
+++ b/test/test_sam.py
@@ -313,7 +313,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testMixedMatch(self):
         """
@@ -328,7 +328,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testHardClipLeft(self):
         """
@@ -343,7 +343,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testHardClipRight(self):
         """
@@ -358,22 +358,23 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testRcNeeded(self):
         """
         A reverse-complimented match (flag = 16) when rcNeeded=True is passed
-        must result in the expected (reverse complimented) padded sequence.
+        must result in the expected (reverse complimented) padded sequence
+        and reversed quality string.
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
-            'query1 16 ref1 2 60 6M * 0 0 TCTAGG ZZZZZZ',
+            'query1 16 ref1 2 60 6M * 0 0 TCTAGG 123456',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries(rcNeeded=True))
-            self.assertEqual(Read('query1', '-CCTAGA---'), read)
+            self.assertEqual(Read('query1', '-CCTAGA---', '!654321!!!'), read)
 
     def testRcSuffix(self):
         """
@@ -382,13 +383,14 @@ class TestPaddedSAM(TestCase):
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
-            'query1 16 ref1 2 60 6M * 0 0 TCTAGG ZZZZZZ',
+            'query1 16 ref1 2 60 6M * 0 0 TCTAGG 123456',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            (read,) = list(ps.queries(rcSuffix='-rc'))
-            self.assertEqual(Read('query1-rc', '-TCTAGG---'), read)
+            (read,) = list(ps.queries(rcSuffix='-rc', rcNeeded=True))
+            self.assertEqual(Read('query1-rc', '-CCTAGA---', '!654321!!!'),
+                             read)
 
     def testQuerySoftClipLeft(self):
         """
@@ -403,7 +405,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testQuerySoftClipReachesLeftEdge(self):
         """
@@ -418,7 +420,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', 'TCTAGG----'), read)
+            self.assertEqual(Read('query1', 'TCTAGG----', 'ZZZZZZ!!!!'), read)
 
     def testQuerySoftClipProtrudesLeft(self):
         """
@@ -433,7 +435,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', 'AGG-------'), read)
+            self.assertEqual(Read('query1', 'AGG-------', 'ZZZ!!!!!!!'), read)
 
     def testKF414679SoftClipLeft(self):
         """
@@ -451,7 +453,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', seq[14:]), read)
+            self.assertEqual(Read('query1', seq[14:], quality[14:]), read)
 
     def testQuerySoftClipRight(self):
         """
@@ -466,7 +468,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '---TCTAGG-'), read)
+            self.assertEqual(Read('query1', '---TCTAGG-', '!!!ZZZZZZ!'), read)
 
     def testQuerySoftClipReachesRightEdge(self):
         """
@@ -481,7 +483,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '----TCTAGG'), read)
+            self.assertEqual(Read('query1', '----TCTAGG', '!!!!ZZZZZZ'), read)
 
     def testQuerySoftClipProtrudesRight(self):
         """
@@ -496,7 +498,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-----TCTAG'), read)
+            self.assertEqual(Read('query1', '-----TCTAG', '!!!!!ZZZZZ'), read)
 
     def testQuerySoftClipProtrudesBothSides(self):
         """
@@ -511,7 +513,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', 'TAGGCTGACT'), read)
+            self.assertEqual(Read('query1', 'TAGGCTGACT', 'ZZZZZZZZZZ'), read)
 
     def testQueryHardClipAndSoftClipProtrudesBothSides(self):
         """
@@ -527,7 +529,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', 'TAGGCTGACT'), read)
+            self.assertEqual(Read('query1', 'TAGGCTGACT', 'ZZZZZZZZZZ'), read)
 
     def testReferenceInsertion(self):
         """
@@ -542,7 +544,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCGG-----'), read)
+            self.assertEqual(Read('query1', '-TCGG-----', '!ZZZZ!!!!!'), read)
             self.assertEqual(
                 {
                     'query1': [(3, 'TA')],
@@ -564,8 +566,9 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read1, read2) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCGG-----'), read1)
-            self.assertEqual(Read('query1/1', '---TCG----'), read2)
+            self.assertEqual(Read('query1', '-TCGG-----', '!ZZZZ!!!!!'), read1)
+            self.assertEqual(Read('query1/1', '---TCG----', '!!!ZZZ!!!!'),
+                             read2)
             self.assertEqual(
                 {
                     'query1': [(3, 'TA')],
@@ -576,7 +579,7 @@ class TestPaddedSAM(TestCase):
     def testReferenceDeletion(self):
         """
         An deletion of reference bases must result in the expected padded
-        sequence (with gaps).
+        sequence (with Ns inserted for the deleted reference bases).
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
@@ -586,12 +589,14 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCNNTAGG-'), read)
+            self.assertEqual(Read('query1', '-TCNNTAGG-', '!ZZ!!ZZZZ!'), read)
 
-    def testReferenceDeletionAlternateChar(self):
+    def testReferenceDeletionAlternateChars(self):
         """
         An deletion of reference bases must result in the expected padded
-        sequence (with gaps) when a queryInsertionChar is passed
+        sequence (with the passed query insertion character and unknown
+        quality character) when queryInsertionChar and unknownQualityChar
+        arguments are passed.
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
@@ -600,13 +605,15 @@ class TestPaddedSAM(TestCase):
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            (read,) = list(ps.queries(queryInsertionChar='?'))
-            self.assertEqual(Read('query1', '-TC??TAGG-'), read)
+            (read,) = list(ps.queries(queryInsertionChar='?',
+                                      unknownQualityChar='+'))
+            self.assertEqual(Read('query1', '-TC??TAGG-', '+ZZ++ZZZZ+'), read)
 
     def testReferenceSkip(self):
         """
         An skip of reference bases must result in the expected padded
-        sequence (with gaps).
+        sequence with the passed unknown quality character when the
+        unknownQualityChar argument is passed.
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
@@ -615,13 +622,15 @@ class TestPaddedSAM(TestCase):
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCNNTAGG-'), read)
+            (read,) = list(ps.queries(unknownQualityChar='.'))
+            self.assertEqual(Read('query1', '-TCNNTAGG-', '.ZZ..ZZZZ.'), read)
 
-    def testReferenceSkipAlternateChar(self):
+    def testReferenceSkipAlternateChars(self):
         """
         An skip of reference bases must result in the expected padded
-        sequence (with gaps) when a queryInsertionChar is passed.
+        sequence (with the passed query insertion character and unknown
+        quality character) when queryInsertionChar and unknownQualityChar
+        arguments are passed.
         """
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
@@ -630,8 +639,9 @@ class TestPaddedSAM(TestCase):
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            (read,) = list(ps.queries(queryInsertionChar='X'))
-            self.assertEqual(Read('query1', '-TCXXTAGG-'), read)
+            (read,) = list(ps.queries(queryInsertionChar='X',
+                                      unknownQualityChar='+'))
+            self.assertEqual(Read('query1', '-TCXXTAGG-', '+ZZ++ZZZZ+'), read)
 
     def testMixedMatchSpecificReferenceButNoMatches(self):
         """
@@ -662,7 +672,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename, referenceIds={'ref1'}))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testMinLength(self):
         """
@@ -679,7 +689,7 @@ class TestPaddedSAM(TestCase):
             filterRead = ReadFilter(minLength=6).filter
             ps = PaddedSAM(SAMFilter(filename, filterRead=filterRead))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testDropSecondary(self):
         """
@@ -694,7 +704,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename, dropSecondary=True))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testDropSupplementary(self):
         """
@@ -710,7 +720,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename, dropSupplementary=True))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testDropDuplicates(self):
         """
@@ -726,7 +736,7 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename, dropDuplicates=True))
             (read,) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read)
 
     def testAllowDuplicateIds(self):
         """
@@ -742,8 +752,9 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read1, read2) = list(ps.queries(allowDuplicateIds=True))
-            self.assertEqual(Read('query1', '-TCTAGG---'), read1)
-            self.assertEqual(Read('query1', '--TC------'), read2)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read1)
+            self.assertEqual(Read('query1', '--TC------', '!!ZZ!!!!!!'),
+                             read2)
 
     def testDuplicateIdDisambiguation(self):
         """
@@ -759,9 +770,11 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read1, read2, read3) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read1)
-            self.assertEqual(Read('query1/1', '--TC------'), read2)
-            self.assertEqual(Read('query1/2', '--TCGA----'), read3)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read1)
+            self.assertEqual(Read('query1/1', '--TC------', '!!ZZ!!!!!!'),
+                             read2)
+            self.assertEqual(Read('query1/2', '--TCGA----', '!!ZZZZ!!!!'),
+                             read3)
 
     def testKeepQualityControlFailures(self):
         """
@@ -777,8 +790,8 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename, keepQCFailures=True))
             (read1, read2) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCTAGG---'), read1)
-            self.assertEqual(Read('query2', '---TC-----'), read2)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read1)
+            self.assertEqual(Read('query2', '---TC-----', '!!!ZZ!!!!!'), read2)
 
     def testSecondaryWithNoPreviousSequence(self):
         """
@@ -792,8 +805,9 @@ class TestPaddedSAM(TestCase):
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            error = ('^Query line 1 has an empty SEQ field, but no previous '
-                     'alignment is present\\.$')
+            error = ('^pysam produced an alignment \\(number 1\\) with no '
+                     'query sequence without previously giving an alignment '
+                     'with a sequence\\.$')
             queries = ps.queries()
             assertRaisesRegex(self, InvalidSAM, error, list, queries)
 
@@ -812,9 +826,10 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read1, read2, read3) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCT------'), read1)
-            self.assertEqual(Read('query2', '-TCTA-----'), read2)
-            self.assertEqual(Read('query2/1', '-----TCTA-'), read3)
+            self.assertEqual(Read('query1', '-TCT------', '!ZZZ!!!!!!'), read1)
+            self.assertEqual(Read('query2', '-TCTA-----', '!ZZZZ!!!!!'), read2)
+            self.assertEqual(Read('query2/1', '-----TCTA-', '!!!!!ZZZZ!'),
+                             read3)
 
     def testSupplementaryWithNoPreviousSequence(self):
         """
@@ -828,8 +843,9 @@ class TestPaddedSAM(TestCase):
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            error = ('^Query line 1 has an empty SEQ field, but no previous '
-                     'alignment is present\\.$')
+            error = ('^pysam produced an alignment \\(number 1\\) with no '
+                     'query sequence without previously giving an alignment '
+                     'with a sequence\\.$')
             queries = ps.queries()
             assertRaisesRegex(self, InvalidSAM, error, list, queries)
 
@@ -848,9 +864,10 @@ class TestPaddedSAM(TestCase):
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read1, read2, read3) = list(ps.queries())
-            self.assertEqual(Read('query1', '-TCT------'), read1)
-            self.assertEqual(Read('query2', '-TCTA-----'), read2)
-            self.assertEqual(Read('query2/1', '-----TCTA-'), read3)
+            self.assertEqual(Read('query1', '-TCT------', '!ZZZ!!!!!!'), read1)
+            self.assertEqual(Read('query2', '-TCTA-----', '!ZZZZ!!!!!'), read2)
+            self.assertEqual(Read('query2/1', '-----TCTA-', '!!!!!ZZZZ!'),
+                             read3)
 
     def testNotSecondaryAndNotSupplementaryWithNoSequence(self):
         """
@@ -864,8 +881,9 @@ class TestPaddedSAM(TestCase):
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            error = ('^Query line 1 has an empty SEQ field, but the alignment '
-                     'is not marked as secondary or supplementary\\.$')
+            error = ('^pysam produced an alignment \\(number 1\\) with no '
+                     'query sequence without previously giving an alignment '
+                     'with a sequence\\.$')
             queries = ps.queries()
             assertRaisesRegex(self, InvalidSAM, error, list, queries)
 
@@ -877,28 +895,28 @@ class TestPaddedSAM(TestCase):
         data = '\n'.join([
             '@SQ SN:ref1 LN:10',
             'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG 123456',
-            'query2 0 ref1 2 60 2= * 0 0 TC XY',
+            'query2 0 ref1 2 60 2= * 0 0 TC 78',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
             (read1, read2) = list(ps.queries(addAlignment=True))
 
-            self.assertEqual(Read('query1', '-TCTAGG---'), read1)
+            self.assertEqual(Read('query1', '-TCTAGG---', '!123456!!!'), read1)
             self.assertEqual('TCTAGG', read1.alignment.query_sequence)
             self.assertEqual('123456', ''.join(
                 map(lambda x: chr(x + 33), read1.alignment.query_qualities)))
 
-            self.assertEqual(Read('query2', '-TC-------'), read2)
+            self.assertEqual(Read('query2', '-TC-------', '!78!!!!!!!'), read2)
             self.assertEqual('TC', read2.alignment.query_sequence)
-            self.assertEqual('XY', ''.join(
+            self.assertEqual('78', ''.join(
                 map(lambda x: chr(x + 33), read2.alignment.query_qualities)))
 
     def testHardClippingInCIGARButQueryNotHardClipped(self):
         """
         As documented in https://github.com/acorg/dark-matter/issues/630 we
-        have to deal correctly with a case in which the CIGAR string says a
-        query should be hard clipped but the query sequence in the SAM file
+        must deal correctly with a case in which the CIGAR string says a
+        query is hard-clipped but the query sequence in the SAM file
         actually isn't. This can be due to a prior alignment with a soft clip,
         in which case the full query sequence has to be given before the
         secondary alignment with the hard clip.
@@ -906,20 +924,77 @@ class TestPaddedSAM(TestCase):
         data = '\n'.join([
             '@SQ SN:Chimp-D00220 LN:8',
             '@SQ SN:D-AM494716 LN:8',
+            '@SQ SN:D-XXX LN:8',
+            '@SQ SN:Chimp-YYY LN:8',
             'query1 0 Chimp-D00220 1 0 3S5M * 0 0 TTTTGGTT 12345678',
             'query1 256 D-AM494716 1 0 3H5M * 0 0 * *',
+            'query1 256 D-XXX 1 0 5H3M * 0 0 * *',
+            'query1 0 Chimp-YYY 1 0 8M * 0 0 * *',
         ]).replace(' ', '\t')
 
         with dataFile(data) as filename:
             ps = PaddedSAM(SAMFilter(filename))
-            (read1, read2) = list(ps.queries(addAlignment=True))
+            (read1, read2, read3, read4) = list(ps.queries(addAlignment=True))
 
-            self.assertEqual(Read('query1', 'TGGTT---'), read1)
+            self.assertEqual(Read('query1', 'TGGTT---', '45678!!!'), read1)
             self.assertEqual('TTTTGGTT', read1.alignment.query_sequence)
 
-            self.assertEqual(Read('query1/1', 'TGGTT---'), read2)
-            # pysam uses None for the query sequence on a secondary alignment.
-            self.assertIs(None, read2.alignment.query_sequence)
+            self.assertEqual(Read('query1/1', 'TGGTT---', '45678!!!'), read2)
+            self.assertEqual('TGGTT', read2.alignment.query_sequence)
+
+            self.assertEqual(Read('query1/2', 'GTT-----', '678!!!!!'), read3)
+            self.assertEqual('GTT', read3.alignment.query_sequence)
+
+            self.assertEqual(Read('query1/3', 'TTTTGGTT', '12345678'), read4)
+            self.assertEqual('TTTTGGTT', read4.alignment.query_sequence)
+
+    def testSecondaryAlignmentHasQuery(self):
+        """
+        If the first alignment of a query is against a reference that is not
+        wanted, a subsequent secondary alignment (SAM flag = 256) must have
+        the original query and quality strings (even though these are only
+        present in the SAM as * characters and the query is None when it comes
+        back from pysam).
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            '@SQ SN:ref2 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query1 256 ref2 2 60 2=2X2M * 0 0 * *',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            ps = PaddedSAM(SAMFilter(filename))
+            (read1, read2) = list(ps.queries(addAlignment=True))
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read1)
+            self.assertEqual('ref1', read1.alignment.reference_name)
+            self.assertEqual(Read('query1/1', '-TCTAGG---', '!ZZZZZZ!!!'),
+                             read2)
+            self.assertEqual('ref2', read2.alignment.reference_name)
+
+    def testSupplementaryAlignmentHasQuery(self):
+        """
+        If the first alignment of a query is against a reference that is not
+        wanted, a subsequent supplementary alignment (SAM flag = 2048) must
+        have the original query and quality strings (even though these are only
+        present in the SAM as * characters and the query is None when it comes
+        back from pysam).
+        """
+        data = '\n'.join([
+            '@SQ SN:ref1 LN:10',
+            '@SQ SN:ref2 LN:10',
+            'query1 0 ref1 2 60 2=2X2M * 0 0 TCTAGG ZZZZZZ',
+            'query1 2048 ref2 2 60 2=2X2M * 0 0 * *',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            ps = PaddedSAM(SAMFilter(filename))
+            (read1, read2) = list(ps.queries(addAlignment=True))
+            self.assertEqual(Read('query1', '-TCTAGG---', '!ZZZZZZ!!!'), read1)
+            self.assertEqual('ref1', read1.alignment.reference_name)
+            self.assertEqual(Read('query1/1', '-TCTAGG---', '!ZZZZZZ!!!'),
+                             read2)
+            self.assertEqual('ref2', read2.alignment.reference_name)
 
 
 class TestSamReferencesToStr(TestCase):
@@ -958,20 +1033,6 @@ class TestHardClip(TestCase):
     """
     Test the _hardClip function.
     """
-    def testCIGARLengthTooHigh(self):
-        """
-        If the total length of the CIGAR operations exceeds the length of the
-        sequence, an AssertionError must be raised.
-        """
-        self.assertRaises(AssertionError, _hardClip, 'CGT', ((CMATCH, 5),))
-
-    def testCIGARLengthTooLow(self):
-        """
-        If the total length of the CIGAR operations is less than the length of
-        the sequence, an AssertionError must be raised.
-        """
-        self.assertRaises(AssertionError, _hardClip, 'CGT', ((CMATCH, 2),))
-
     def testHardClipInMiddle(self):
         """
         If hard clipping is given as an operation not at the beginning or end
@@ -980,8 +1041,8 @@ class TestHardClip(TestCase):
         error = ('^Invalid CIGAR tuples .* contains hard-clipping operation '
                  'that is neither at the start nor the end of the sequence\.$')
         self.assertRaisesRegex(
-            ValueError, error,
-            _hardClip, 'CGT', ((CMATCH, 1), (CHARD_CLIP, 1), (CMATCH, 1),))
+            ValueError, error, _hardClip, 'CGT', '123',
+            ((CMATCH, 1), (CHARD_CLIP, 1), (CMATCH, 1),))
 
     def testThreeHardClips(self):
         """
@@ -991,8 +1052,7 @@ class TestHardClip(TestCase):
         error = ('^Invalid CIGAR tuples .* specifies hard-clipping 3 times '
                  '\(2 is the maximum\).$')
         self.assertRaisesRegex(
-            ValueError, error,
-            _hardClip, 'CGT',
+            ValueError, error, _hardClip, 'CGT', '123',
             ((CHARD_CLIP, 1), (CHARD_CLIP, 1), (CHARD_CLIP, 1),))
 
     def testNoClip(self):
@@ -1000,23 +1060,26 @@ class TestHardClip(TestCase):
         If no hard clipping is indicated, the function must return the
         original sequence.
         """
-        self.assertEqual('CGT', _hardClip('CGT', ((CMATCH, 3),)))
+        self.assertEqual(('CGT', '123', False),
+                         _hardClip('CGT', '123', ((CMATCH, 3),)))
 
     def testClipLeft(self):
         """
         If hard clipping on the left is indicated, and has not been done,
         the function must return the expected sequence.
         """
-        self.assertEqual('CGT',
-                         _hardClip('CAACGT', ((CHARD_CLIP, 3), (CMATCH, 3),)))
+        self.assertEqual(
+            ('CGT', '456', True),
+            _hardClip('CAACGT', '123456', ((CHARD_CLIP, 3), (CMATCH, 3),)))
 
     def testClipRight(self):
         """
         If hard clipping on the right is indicated, and has not been done,
         the function must return the expected sequence.
         """
-        self.assertEqual('CA',
-                         _hardClip('CAACGT', ((CMATCH, 2), (CHARD_CLIP, 4),)))
+        self.assertEqual(
+            ('CA', '12', True),
+            _hardClip('CAACGT', '123456', ((CMATCH, 2), (CHARD_CLIP, 4),)))
 
     def testClipBoth(self):
         """
@@ -1024,8 +1087,8 @@ class TestHardClip(TestCase):
         done, the function must return the expected sequence.
         """
         self.assertEqual(
-            'AA',
-            _hardClip('CAACGT',
+            ('AA', '23', True),
+            _hardClip('CAACGT', '123456',
                       ((CHARD_CLIP, 1), (CMATCH, 2), (CHARD_CLIP, 3),)))
 
     def testClipLeftAlreadyDone(self):
@@ -1033,16 +1096,18 @@ class TestHardClip(TestCase):
         If hard clipping on the left is indicated, and has already been done,
         the function must return the expected sequence.
         """
-        self.assertEqual('CGT',
-                         _hardClip('CGT', ((CHARD_CLIP, 3), (CMATCH, 3),)))
+        self.assertEqual(
+            ('CGT', '123', False),
+            _hardClip('CGT', '123', ((CHARD_CLIP, 3), (CMATCH, 3),)))
 
     def testClipRightAlreadyDone(self):
         """
         If hard clipping on the right is indicated, and has already been done,
         the function must return the expected sequence.
         """
-        self.assertEqual('CA',
-                         _hardClip('CA', ((CMATCH, 2), (CHARD_CLIP, 4),)))
+        self.assertEqual(
+            ('CA', '12', False),
+            _hardClip('CA', '12', ((CMATCH, 2), (CHARD_CLIP, 4),)))
 
     def testClipBothAlreadyDone(self):
         """
@@ -1050,6 +1115,6 @@ class TestHardClip(TestCase):
         been done, the function must return the expected sequence.
         """
         self.assertEqual(
-            'AA',
-            _hardClip('AA',
+            ('AA', '12', False),
+            _hardClip('AA', '12',
                       ((CHARD_CLIP, 1), (CMATCH, 2), (CHARD_CLIP, 3),)))