Repository: acorg/dark-matter

Issue:
SAMFilter raises if pysam returns an alignment with a None quality string
If a SAM file has `*` for a quality string, `pysam` gives back a `None` field and this causes an exception in `dark/sam.py`:

```
$ read-sam.py 180426-18-mouse-L1.sam
Traceback (most recent call last):
  File "./read-sam.py", line 8, in <module>
    for alignment in sam.alignments():
  File "/Users/terry/dark-matter/dark-matter/dark/sam.py", line 365, in alignments
    alignment.cigartuples)
  File "/Users/terry/dark-matter/dark-matter/dark/sam.py", line 93, in _hardClip
    assert sequenceLength == len(quality)
TypeError: object of type 'NoneType' has no len()
```

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 588c254..8ddbc7f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,8 @@
+## 3.0.53 Nov 23, 2018
+
+Fixed [#650](https://github.com/acorg/dark-matter/issues/650),
+exception in `SAMFilter` when quality is `*` in a SAM file.
+
 ## 3.0.52 Nov 23, 2018
 
 Added hard-clipping to CIGAR in SAM created by `convert-diamond-to-sam.py`.
diff --git a/dark/__init__.py b/dark/__init__.py
index 147b204..2f33ce1 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.52'
+__version__ = '3.0.53'
diff --git a/dark/sam.py b/dark/sam.py
index f6052ca..b86dbb2 100644
--- a/dark/sam.py
+++ b/dark/sam.py
@@ -70,8 +70,9 @@ def _hardClip(sequence, quality, cigartuples):
     Hard clip (if necessary) a sequence.
 
     @param sequence: A C{str} nucleotide sequence.
-    @param quality: A C{str} quality string or a C{list} of C{int} quality
-        values as returned by pysam.
+    @param quality: A C{str} quality string, or a C{list} of C{int} quality
+        values as returned by pysam, or C{None} if the SAM file had a '*'
+        for the quality string (which pysam converts to C{None}).
     @param cigartuples: An iterable of (operation, length) tuples, detailing
         the alignment, as per the SAM specification.
     @return: A 3-tuple consisting of
@@ -90,7 +91,8 @@ def _hardClip(sequence, quality, cigartuples):
         cigarLength += length if operation in _CONSUMES_QUERY else 0
 
     sequenceLength = len(sequence)
-    assert sequenceLength == len(quality)
+    if quality is not None:
+        assert sequenceLength == len(quality)
     clipLeft = clipRight = 0
     clippedSequence = sequence
     clippedQuality = quality
@@ -109,12 +111,14 @@ def _hardClip(sequence, quality, cigartuples):
             if not alreadyClipped:
                 clipLeft = cigartuples[0][1]
                 clippedSequence = sequence[clipLeft:]
-                clippedQuality = quality[clipLeft:]
+                if quality is not None:
+                    clippedQuality = quality[clipLeft:]
         elif cigartuples[-1][0] == CHARD_CLIP:
             if not alreadyClipped:
                 clipRight = cigartuples[-1][1]
                 clippedSequence = sequence[:-clipRight]
-                clippedQuality = quality[:-clipRight]
+                if quality is not None:
+                    clippedQuality = quality[:-clipRight]
         else:
             raise ValueError(
                 'Invalid CIGAR tuples (%s) contains hard-clipping operation '
@@ -126,7 +130,8 @@ def _hardClip(sequence, quality, cigartuples):
         if not alreadyClipped:
             clipLeft, clipRight = cigartuples[0][1], cigartuples[-1][1]
             clippedSequence = sequence[clipLeft:-clipRight]
-            clippedQuality = quality[clipLeft:-clipRight]
+            if quality is not None:
+                clippedQuality = quality[clipLeft:-clipRight]
     else:
         raise ValueError(
             'Invalid CIGAR tuples (%s) specifies hard-clipping %d times (2 '
@@ -147,7 +152,9 @@ def _hardClip(sequence, quality, cigartuples):
                  abs(len(sequence) - len(clippedSequence)),
                  clipLeft + clipRight, clipLeft, clipRight, cigartuples))
     else:
-        assert len(clippedSequence) == len(clippedQuality) == sequenceLength
+        assert len(clippedSequence) == sequenceLength
+        if quality is not None:
+            assert len(clippedQuality) == sequenceLength
 
     return clippedSequence, clippedQuality, weClipped
 


Test patch:
diff --git a/test/test_sam.py b/test/test_sam.py
index ca2e4d0..42ef202 100644
--- a/test/test_sam.py
+++ b/test/test_sam.py
@@ -265,6 +265,23 @@ class TestSAMFilter(TestCase):
             self.assertEqual('ZZZZZZ', ''.join(
                 map(lambda x: chr(x + 33), alignment.query_qualities)))
 
+    def testNoQuality(self):
+        """
+        If an alignment has * for the quality string, the filter must
+        return an alignment with a C{None} quality value.
+        """
+        data = '\n'.join([
+            '@SQ SN:ref LN:10',
+            'query1 4 * 0 0 6M * 0 0 TCTAGG *',
+        ]).replace(' ', '\t')
+
+        with dataFile(data) as filename:
+            sf = SAMFilter(filename)
+            (alignment,) = list(sf.alignments())
+            self.assertEqual('query1', alignment.query_name)
+            self.assertEqual('TCTAGG', alignment.query_sequence)
+            self.assertIsNone(alignment.query_qualities)
+
 
 class TestPaddedSAM(TestCase):
     """