Repository: acorg/dark-matter

Issue:
Add soft-clipping to CIGAR string in convert-diamond-to-sam.py

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 271f053..588c254 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,7 @@
+## 3.0.52 Nov 23, 2018
+
+Added hard-clipping to CIGAR in SAM created by `convert-diamond-to-sam.py`.
+
 ## 3.0.51 Nov 23, 2018
 
 Use `from six import StringIO` to avoid a PY2/3 incompatibility.
diff --git a/bin/convert-diamond-to-sam.py b/bin/convert-diamond-to-sam.py
index 38e3365..0b3484f 100755
--- a/bin/convert-diamond-to-sam.py
+++ b/bin/convert-diamond-to-sam.py
@@ -16,7 +16,8 @@ from dark.btop import btop2cigar
 from dark.diamond.conversion import diamondTabularFormatToDicts
 from dark.reads import DNARead
 
-FIELDS = 'bitscore btop qframe qqual qseq qseqid slen sstart stitle'
+FIELDS = ('bitscore btop qframe qend qqual qlen qseq qseqid qstart slen '
+          'sstart stitle')
 
 parser = argparse.ArgumentParser(
     formatter_class=argparse.ArgumentDefaultsHelpFormatter,
@@ -92,6 +93,21 @@ for match in diamondTabularFormatToDicts(sys.stdin, FIELDS.split()):
         qseq = DNARead('id', match['qseq']).reverseComplement().sequence
         qqual = match['qqual'][::-1] if match['qqual'] else '*'
 
+    # Make a CIGAR string, including hard-clipped bases at the start and
+    # end of the query (DIAMOND outputs a hard-clipped query sequence).
+    startClipCount = match['qstart'] - 1
+    endClipCount = match['qlen'] - match['qend']
+
+    assert startClipCount >= 0
+    assert endClipCount >= 0, (
+        'Query sequence %s has length %d but the qend value is %d' %
+        (qseq, len(match['qseq']), match['qend']))
+
+    cigar = (
+        ('%dH' % startClipCount if startClipCount else '') +
+        btop2cigar(match['btop'], concise=False, aa=True) +
+        ('%dH' % endClipCount if endClipCount else ''))
+
     emit('\t'.join(map(str, [
         # 1. QNAME
         qseqid,
@@ -104,7 +120,7 @@ for match in diamondTabularFormatToDicts(sys.stdin, FIELDS.split()):
         # 5. MAPQ
         mappingQuality,
         # 6. CIGAR
-        btop2cigar(match['btop']),
+        cigar,
         # 7. RNEXT
         '*',
         # 8. PNEXT
diff --git a/dark/__init__.py b/dark/__init__.py
index 3e72b1a..147b204 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.51'
+__version__ = '3.0.52'
diff --git a/dark/btop.py b/dark/btop.py
index 9855734..aa338da 100644
--- a/dark/btop.py
+++ b/dark/btop.py
@@ -73,24 +73,34 @@ def countGaps(btopString):
     return (queryGaps, subjectGaps)
 
 
-def btop2cigar(btopString, concise=False):
+def btop2cigar(btopString, concise=False, aa=False):
     """
     Convert a BTOP string to a CIGAR string.
 
     @param btopString: A C{str} BTOP sequence.
     @param concise: If C{True}, use 'M' for matches and mismatches instead
         of the more specific 'X' and '='.
-
-    @raise ValueError: If L{parseBtop} finds an error in C{btopString}.
+    @param aa: If C{True}, C{btopString} will be interpreted as though it
+        refers to amino acids (as in the BTOP string produced by DIAMOND).
+        In that case, it is not possible to use the 'precise' CIGAR characters
+        because amino acids have multiple codons so we cannot know whether
+        an amino acid match is due to an exact nucleotide matches or not.
+        Also, the numbers in the BTOP string will be multiplied by 3 since
+        they refer to a number of amino acids matching.
+    @raise ValueError: If L{parseBtop} finds an error in C{btopString} or
+        if C{aa} and C{concise} are both C{True}.
     @return: A C{str} CIGAR string.
     """
+    if aa and concise:
+        raise ValueError('aa and concise cannot both be True')
+
     result = []
     thisLength = thisOperation = currentLength = currentOperation = None
 
     for item in parseBtop(btopString):
         if isinstance(item, int):
             thisLength = item
-            thisOperation = CMATCH if concise else CEQUAL
+            thisOperation = CEQUAL if concise else CMATCH
         else:
             thisLength = 1
             query, reference = item
@@ -108,13 +118,16 @@ def btop2cigar(btopString, concise=False):
             else:
                 # A substitution was needed.
                 assert query != reference
-                thisOperation = CMATCH if concise else CDIFF
+                thisOperation = CDIFF if concise else CMATCH
 
         if thisOperation == currentOperation:
             currentLength += thisLength
         else:
             if currentOperation:
-                result.append('%d%s' % (currentLength, currentOperation))
+                result.append(
+                    '%d%s' %
+                    ((3 * currentLength) if aa else currentLength,
+                     currentOperation))
             currentLength, currentOperation = thisLength, thisOperation
 
     # We reached the end of the BTOP string. If there was an operation
@@ -122,6 +135,8 @@ def btop2cigar(btopString, concise=False):
     # case where btopString was empty.
     assert currentOperation or btopString == ''
     if currentOperation:
-        result.append('%d%s' % (currentLength, currentOperation))
+        result.append(
+            '%d%s' %
+            ((3 * currentLength) if aa else currentLength, currentOperation))
 
     return ''.join(result)


Test patch:
diff --git a/test/test_btop.py b/test/test_btop.py
index 9f1f22a..e019563 100644
--- a/test/test_btop.py
+++ b/test/test_btop.py
@@ -164,15 +164,15 @@ class TestCountGaps(TestCase):
         self.assertEqual((3, 2), countGaps('-GG-34-T-T39F-'))
 
 
-class TestBtop2CigarConcise(TestCase):
+class TestBtop2CigarPrecise(TestCase):
     """
-    Tests for the btop2cigar function when concise is True.
+    Tests for the btop2cigar function when concise is False.
     """
     def testEmpty(self):
         """
         An empty BTOP string must result in an empty CIGAR string.
         """
-        self.assertEqual('', btop2cigar('', concise=True))
+        self.assertEqual('', btop2cigar('', concise=False))
 
     def testMixedMatch(self):
         """
@@ -180,7 +180,7 @@ class TestBtop2CigarConcise(TestCase):
         CIGAR sense where M could be identical characters or not), the CIGAR
         string must be all Ms.
         """
-        self.assertEqual('7M', btop2cigar('2GC3AT', concise=True))
+        self.assertEqual('7M', btop2cigar('2GC3AT', concise=False))
 
     def testRefenceInsertion(self):
         """
@@ -188,7 +188,7 @@ class TestBtop2CigarConcise(TestCase):
         subject (reference) does not, the CIGAR string must indicate an
         insertion to the reference.
         """
-        self.assertEqual('1I', btop2cigar('A-', concise=True))
+        self.assertEqual('1I', btop2cigar('A-', concise=False))
 
     def testQueryInsertion(self):
         """
@@ -196,7 +196,7 @@ class TestBtop2CigarConcise(TestCase):
         but the query does not, the CIGAR string must indicate an deletion in
         the reference.
         """
-        self.assertEqual('1D', btop2cigar('-A', concise=True))
+        self.assertEqual('1D', btop2cigar('-A', concise=False))
 
     def testAll(self):
         """
@@ -204,25 +204,36 @@ class TestBtop2CigarConcise(TestCase):
         string must be correct.
         """
         self.assertEqual('7M2I4M2D5M',
-                         btop2cigar('2GC3ATC-G-4-T-A5', concise=True))
+                         btop2cigar('2GC3ATC-G-4-T-A5', concise=False))
 
+    def testAllAA(self):
+        """
+        If a BTOP string specifies all possible variations, and we indicate
+        that the BTOP string refers to amino acids, the CIGAR string must be
+        correct (i.e., all counts must be tripled).
+        """
+        self.assertEqual(
+            '21M6I12M6D15M',
+            btop2cigar('2GC3ATC-G-4-T-A5', concise=False, aa=True))
 
-class TestBtop2CigarPrecise(TestCase):
+
+class TestBtop2CigarConcise(TestCase):
     """
-    Tests for the btop2cigar function when concise is False.
+    Tests for the btop2cigar function when concise is True.
     """
+
     def testEmpty(self):
         """
         An empty BTOP string must result in an empty CIGAR string.
         """
-        self.assertEqual('', btop2cigar('', concise=False))
+        self.assertEqual('', btop2cigar('', concise=True))
 
     def testMixedMatch(self):
         """
         If a BTOP string specifies that some characters match and some do
         not, the CIGAR string must be specific about the matches / mismatches.
         """
-        self.assertEqual('2=1X3=1X', btop2cigar('2GC3AT', concise=False))
+        self.assertEqual('2=1X3=1X', btop2cigar('2GC3AT', concise=True))
 
     def testRefenceInsertion(self):
         """
@@ -230,7 +241,7 @@ class TestBtop2CigarPrecise(TestCase):
         subject (reference) does not, the CIGAR string must indicate an
         insertion to the reference.
         """
-        self.assertEqual('1I', btop2cigar('A-', concise=False))
+        self.assertEqual('1I', btop2cigar('A-', concise=True))
 
     def testQueryInsertion(self):
         """
@@ -238,7 +249,7 @@ class TestBtop2CigarPrecise(TestCase):
         but the query does not, the CIGAR string must indicate an deletion in
         the reference.
         """
-        self.assertEqual('1D', btop2cigar('-A', concise=False))
+        self.assertEqual('1D', btop2cigar('-A', concise=True))
 
     def testAll(self):
         """
@@ -246,4 +257,12 @@ class TestBtop2CigarPrecise(TestCase):
         string must be correct.
         """
         self.assertEqual('2=1X3=1X2I4=2D5=',
-                         btop2cigar('2GC3ATC-G-4-T-A5', concise=False))
+                         btop2cigar('2GC3ATC-G-4-T-A5', concise=True))
+
+    def testWithAATrue(self):
+        """
+        If concise and aa are both set to True, a ValueError must be raised.
+        """
+        error = '^aa and concise cannot both be True$'
+        assertRaisesRegex(self, ValueError, error, btop2cigar, '',
+                          concise=True, aa=True)