Repository: acorg/dark-matter

Issue:
Add --reverse and --reverseCompliment options to filter-fasta.py

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e329501..2c1e45d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,8 @@
+## 3.0.61 Dec 29, 2018
+
+Added `--reverse` and `--reverseComplement` options to `filter-fasta.py`
+and the underlying `ReadFilter` class.
+
 ## 3.0.60 Dec 13, 2018
 
 In `reads.py`, changed the `_makeComplementTable` function so that 
diff --git a/bin/convert-diamond-to-sam.py b/bin/convert-diamond-to-sam.py
index a4a7bd5..fbab91f 100755
--- a/bin/convert-diamond-to-sam.py
+++ b/bin/convert-diamond-to-sam.py
@@ -77,7 +77,7 @@ for match in diamondTabularFormatToDicts(sys.stdin):
     referenceLengths[stitle] = 3 * match['slen']
 
     # If the query frame is less than zero, the match was with a reverse
-    # complemented translation of the query. Put the reverse compliment
+    # complemented translation of the query. Put the reverse complement
     # into the SAM output, which seems to be standard / accepted practice
     # based on my web searches. See e.g., https://www.biostars.org/p/131891/
     # for what Bowtie2 does and for some comments on this issue for SAM/BAM
diff --git a/dark/__init__.py b/dark/__init__.py
index ffb465c..7f64217 100644
--- a/dark/__init__.py
+++ b/dark/__init__.py
@@ -8,4 +8,4 @@ if sys.version_info < (2, 7):
 #
 # Remember to update ../CHANGELOG.md describing what's new in each version.
 
-__version__ = '3.0.60'
+__version__ = '3.0.61'
diff --git a/dark/filter.py b/dark/filter.py
index 05d36e8..91971bd 100644
--- a/dark/filter.py
+++ b/dark/filter.py
@@ -379,6 +379,15 @@ def addFASTAEditingCommandLineOptions(parser):
               'function will be passed to eval with the dark.reads classes '
               'Read, DNARead, AARead, etc. all in scope.'))
 
+    parser.add_argument(
+        '--reverse', action='store_true', default=False,
+        help=('Reverse the sequences. Note that this is NOT reverse '
+              'complementing.'))
+
+    parser.add_argument(
+        '--reverseComplement', action='store_true', default=False,
+        help='Reverse complement the sequences.')
+
 
 def parseFASTAEditingCommandLineOptions(args, reads):
     """
@@ -430,4 +439,5 @@ def parseFASTAEditingCommandLineOptions(args, reads):
         truncateTitlesAfter=truncateTitlesAfter,
         removeDescriptions=removeDescriptions,
         idLambda=args.idLambda, readLambda=args.readLambda,
-        keepSites=keepSites, removeSites=removeSites)
+        keepSites=keepSites, removeSites=removeSites,
+        reverse=args.reverse, reverseComplement=args.reverseComplement)
diff --git a/dark/reads.py b/dark/reads.py
index df381e0..98a78cc 100644
--- a/dark/reads.py
+++ b/dark/reads.py
@@ -126,6 +126,17 @@ class Read(object):
             'quality': self.quality,
         }
 
+    def reverse(self):
+        """
+        Reverse a read (note that this is NOT a reverse complement).
+
+        @return: The reversed sequence as an instance of the current class.
+        """
+        return self.__class__(
+            self.id,
+            self.sequence[::-1],
+            None if self.quality is None else self.quality[::-1])
+
     @classmethod
     def fromDict(cls, d):
         """
@@ -807,12 +818,20 @@ class ReadFilter(object):
     @param removeSites: A set of C{int} 0-based sites (i.e., indices) in
         sequences that should be removed. If C{None} (the default), no sites
         are removed.
+    @param reverse: If C{True}, reverse the sequences. Reversing happens
+        at a very late stage (i.e., after sites are altered via keepSites
+        and removeSites).
+    @param reverseComplement: If C{True}, replace seqeunces with their reverse
+        complements. Reversing happens at a very late stage (i.e., after sites
+        are altered via keepSites and removeSites).
     @raises ValueError: If C{randomSubset} and C{sampleFraction} are both
         specified, or if C{randomSubset} is specified but C{trueLength} is not,
         or if the sequence numbers in C{sequenceNumbersFile} are
         non-positive or not ascending, or if both C{keepSites} and
         C{removeSites} are given, or if both C{keepSequences} and
         C{removeSequences} are given.
+    @raise AttributeError: If C{reverseComplement} is C{True} but the read
+        type does not allow for reverse complementing.
     """
 
     # TODO, when/if needed: make it possible to pass a seed for the RNG
@@ -830,7 +849,8 @@ class ReadFilter(object):
                  removeDescriptions=False, modifier=None, randomSubset=None,
                  trueLength=None, sampleFraction=None,
                  sequenceNumbersFile=None, idLambda=None, readLambda=None,
-                 keepSites=None, removeSites=None):
+                 keepSites=None, removeSites=None, reverse=False,
+                 reverseComplement=False):
 
         if randomSubset is not None:
             if sampleFraction is not None:
@@ -867,6 +887,12 @@ class ReadFilter(object):
         self.keepSites = keepSites
         self.removeSites = removeSites
 
+        if reverseComplement:
+            # Make sure reverse is not also set.
+            reverse = False
+        self.reverse = reverse
+        self.reverseComplement = reverseComplement
+
         self.alwaysFalse = False
         self.yieldCount = 0
         self.readIndex = -1
@@ -1069,6 +1095,11 @@ class ReadFilter(object):
         if self.removeDescriptions:
             read.id = read.id.split()[0]
 
+        if self.reverse:
+            read = read.reverse()
+        elif self.reverseComplement:
+            read = read.reverseComplement()
+
         self.yieldCount += 1
         return read
 
diff --git a/dark/sequence.py b/dark/sequence.py
index 6266569..028bd00 100644
--- a/dark/sequence.py
+++ b/dark/sequence.py
@@ -30,7 +30,7 @@ def findPrimer(primer, seq):
 
 def findPrimerBidi(primer, seq):
     """
-    Look for a primer in a sequence and its reverse compliment.
+    Look for a primer in a sequence and its reverse complement.
 
     @param primer: A C{str} primer sequence.
     @param seq: A BioPython C{Bio.Seq} sequence.
@@ -38,7 +38,7 @@ def findPrimerBidi(primer, seq):
     @return: A C{tuple} of two lists. The first contains (zero-based)
         ascending offsets into the sequence at which the primer can be found.
         The second is a similar list ascending offsets into the original
-        sequence where the primer matches the reverse complimented of the
+        sequence where the primer matches the reverse complemented of the
         sequence. If no instances are found, the corresponding list in the
         returned tuple must be empty.
     """
@@ -54,7 +54,7 @@ def findPrimerBidi(primer, seq):
 def findPrimerBidiLimits(primer, seq):
     """
     Report the extreme (inner) offsets of primer in a sequence and its
-    reverse compliment.
+    reverse complement.
 
     @param primer: A C{str} primer sequence.
     @param seq: A BioPython C{Bio.Seq} sequence.
diff --git a/requirements.txt b/requirements.txt
index 1c46c1e..8ae39e2 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,7 +10,7 @@ numpy>=1.14.2
 pyfaidx==0.5.0
 pysam>=0.14.1
 pysam>=0.14.1
-pytest>=3.8.1
+pytest>=4.0.2
 pyzmq>=14.3.1
 requests>=2.18.4
 simplejson>=3.5.3


Test patch:
diff --git a/test/test_reads.py b/test/test_reads.py
index 4663234..8339e8a 100644
--- a/test/test_reads.py
+++ b/test/test_reads.py
@@ -619,6 +619,22 @@ class TestRead(TestCase):
                          Read('id1', 'ATCGAT').newFromSites({100, 200},
                                                             exclude=True))
 
+    def testReverseNoQuality(self):
+        """
+        The reverse method must work as expected when there is no quality
+        string.
+        """
+        self.assertEqual(Read('id1', 'ATCGAT'),
+                         Read('id1', 'TAGCTA').reverse())
+
+    def testReverseWithQuality(self):
+        """
+        The reverse method must work as expected when there is a quality
+        string.
+        """
+        self.assertEqual(Read('id1', 'ATCGAT', '123456'),
+                         Read('id1', 'TAGCTA', '654321').reverse())
+
 
 class TestDNARead(TestCase):
     """
@@ -3184,6 +3200,63 @@ class TestReadsFiltering(TestCase):
         (result,) = list(result)
         self.assertEqual(Read('xid1-x', 'AT'), result)
 
+    def testReverse(self):
+        """
+        When reverse=True, reads with the expected sequences and qualities
+        must be returned.
+        """
+        read1 = Read('id1', 'ATCGCC', '123456')
+        read2 = Read('id2', 'GGATCG', '987654')
+        reads = Reads(initialReads=[read1, read2])
+        result1, result2 = list(reads.filter(reverse=True))
+        self.assertEqual(Read('id1', 'CCGCTA', '654321'), result1)
+        self.assertEqual(Read('id2', 'GCTAGG', '456789'), result2)
+
+    def testReverseComplement(self):
+        """
+        When reverseComplement=True, reads with the expected sequences and
+        qualities must be returned when the reads are of type DNARead.
+        """
+        read1 = DNARead('id1', 'ATCGCC', '123456')
+        read2 = DNARead('id2', 'GGATCG', '987654')
+        reads = Reads(initialReads=[read1, read2])
+        result1, result2 = list(reads.filter(reverseComplement=True))
+        self.assertEqual(Read('id1', 'GGCGAT', '654321'), result1)
+        self.assertEqual(Read('id2', 'CGATCC', '456789'), result2)
+
+    def testReverseAndReverseComplement(self):
+        """
+        When reverseComplement=True and reverse=True, the expected reads must
+        be returned (reverse complemented) when the reads are of type DNARead.
+        """
+        read1 = DNARead('id1', 'ATCGCC', '123456')
+        read2 = DNARead('id2', 'GGATCG', '987654')
+        reads = Reads(initialReads=[read1, read2])
+        result1, result2 = list(reads.filter(reverseComplement=True,
+                                             reverse=True))
+        self.assertEqual(Read('id1', 'GGCGAT', '654321'), result1)
+        self.assertEqual(Read('id2', 'CGATCC', '456789'), result2)
+
+    def testReverseComplementNonDNA(self):
+        """
+        When reverseComplement=True and the read is not a DNARead, an
+        AttributeError must be raised.
+        """
+        reads = Reads(initialReads=[Read('id1', 'ATCGCC', '123456')])
+        error = "^'Read' object has no attribute 'reverseComplement'$"
+        six.assertRaisesRegex(self, AttributeError, error, list,
+                              reads.filter(reverseComplement=True))
+
+    def testReverseComplementAARead(self):
+        """
+        When reverseComplement=True and the read is an AARead, an
+        AttributeError must be raised.
+        """
+        reads = Reads(initialReads=[AARead('id1', 'ATCGCC', '123456')])
+        error = "^'AARead' object has no attribute 'reverseComplement'$"
+        six.assertRaisesRegex(self, AttributeError, error, list,
+                              reads.filter(reverseComplement=True))
+
 
 class TestReadsInRAM(TestCase):
     """
diff --git a/test/test_sam.py b/test/test_sam.py
index 42ef202..dc5908a 100644
--- a/test/test_sam.py
+++ b/test/test_sam.py
@@ -397,8 +397,8 @@ class TestPaddedSAM(TestCase):
 
     def testRcNeeded(self):
         """
-        A reverse-complimented match (flag = 16) when rcNeeded=True is passed
-        must result in the expected (reverse complimented) padded sequence
+        A reverse-complemented match (flag = 16) when rcNeeded=True is passed
+        must result in the expected (reverse complemented) padded sequence
         and reversed quality string.
         """
         data = '\n'.join([
@@ -413,7 +413,7 @@ class TestPaddedSAM(TestCase):
 
     def testRcSuffix(self):
         """
-        A reverse-complimented sequence should have the rcSuffix string added
+        A reverse-complemented sequence should have the rcSuffix string added
         to its id when an rcSuffix value is passed.
         """
         data = '\n'.join([