Repository: acorg/dark-matter

Issue:
Add filtering options for maximum number of Ns in a sequence

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 66b60e7..17a90c5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,7 @@
+## 3.1.70 March 24, 2020
+
+Added `--maxNFraction` arg to `filter-fasta.py`.
+
 ## 3.1.69 March 23, 2020
 
 Added `--md5` arg to `fasta-sequences.py`.
@@ -429,9 +433,9 @@ Removed bone-headed use of full path to `fasta-join.sh` from
 
 ## 3.0.63 Jan 14, 2019
 
-Added `compareAaReads` and `matchToString` to `aa.py`. Wrote tests in 
+Added `compareAaReads` and `matchToString` to `aa.py`. Wrote tests in
 `test_aa.py` for both. Moved `countPrint` to utils, used by `matchToString`
-in `dna.py` and `aa.py`. Added `compare-aa-sequences.py` to the bin. 
+in `dna.py` and `aa.py`. Added `compare-aa-sequences.py` to the bin.
 
 ## 3.0.62 Dec 30, 2018
 
@@ -444,9 +448,9 @@ and the underlying `ReadFilter` class.
 
 ## 3.0.60 Dec 13, 2018
 
-In `reads.py`, changed the `_makeComplementTable` function so that 
+In `reads.py`, changed the `_makeComplementTable` function so that
 uppercase and lowercase bases are correctly reverse complemented into their
-respective uppercase and lowercase complementary letters. Added a test to 
+respective uppercase and lowercase complementary letters. Added a test to
 `test/reads.py` to confirm that `reverseComplement` does this.
 
 ## 3.0.59 Dec 11, 2018
diff --git a/dark/__init__.py b/dark/__init__.py
index 52bc094..01d72bb 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.1.69'
+__version__ = '3.1.70'
diff --git a/dark/filter.py b/dark/filter.py
index 07e00d7..216db82 100644
--- a/dark/filter.py
+++ b/dark/filter.py
@@ -1,3 +1,5 @@
+from __future__ import division
+
 import re
 from math import ceil
 from collections import OrderedDict
@@ -187,6 +189,10 @@ def addFASTAFilteringCommandLineOptions(parser):
         '--maxLength', type=int, metavar='N',
         help='The maximum sequence length')
 
+    parser.add_argument(
+        '--maxNFraction', type=float, metavar='N',
+        help='The maximum fraction of Ns that can be present in the sequence')
+
     parser.add_argument(
         '--whitelist', action='append', metavar='SEQUENCE-ID',
         help='Sequence titles (ids) that should be whitelisted')
@@ -297,6 +303,7 @@ def parseFASTAFilteringCommandLineOptions(args, reads):
 
     return reads.filter(
         minLength=args.minLength, maxLength=args.maxLength,
+        maxNFraction=args.maxNFraction,
         whitelist=set(args.whitelist) if args.whitelist else None,
         blacklist=set(args.blacklist) if args.blacklist else None,
         whitelistFile=args.whitelistFile, blacklistFile=args.blacklistFile,
diff --git a/dark/reads.py b/dark/reads.py
index 8da2fef..113ffdd 100644
--- a/dark/reads.py
+++ b/dark/reads.py
@@ -817,6 +817,8 @@ class ReadFilter(object):
 
     @param minLength: The minimum acceptable length.
     @param maxLength: The maximum acceptable length.
+    @param maxNFraction: The maximum fraction of Ns that can be present in the
+        sequence.
     @param removeGaps: If C{True} remove all gaps ('-' characters) from the
         read sequences.
     @param whitelist: If not C{None}, a set of exact read ids that are
@@ -935,8 +937,8 @@ class ReadFilter(object):
     # save and restore the state of the RNG and/or to optionally add
     # 'seed=XXX' to the end of the id of the first read, etc.
 
-    def __init__(self, minLength=None, maxLength=None, removeGaps=False,
-                 whitelist=None, blacklist=None,
+    def __init__(self, minLength=None, maxLength=None, maxNFraction=None,
+                 removeGaps=False, whitelist=None, blacklist=None,
                  whitelistFile=None, blacklistFile=None,
                  titleRegex=None, negativeTitleRegex=None,
                  truncateTitlesAfter=None, keepSequences=None,
@@ -960,6 +962,7 @@ class ReadFilter(object):
 
         self.minLength = minLength
         self.maxLength = maxLength
+        self.maxNFraction = maxNFraction
         self.removeGaps = removeGaps
         self.head = head
         self.removeDuplicates = removeDuplicates
@@ -1126,6 +1129,11 @@ class ReadFilter(object):
                 (self.maxLength is not None and readLen > self.maxLength)):
             return False
 
+        if self.maxNFraction is not None:
+            nFraction = read.sequence.count('N') / readLen
+            if self.maxNFraction < nFraction:
+                return False
+
         if self.removeGaps:
             if read.quality is None:
                 read = read.__class__(read.id, read.sequence.replace('-', ''))


Test patch:
diff --git a/test/test_reads.py b/test/test_reads.py
index 273b87d..ea53b3f 100644
--- a/test/test_reads.py
+++ b/test/test_reads.py
@@ -2550,6 +2550,46 @@ class TestReads(TestCase):
         self.assertEqual(sorted([read2, read3]), sorted(reads))
         self.assertEqual(3, reads.unfilteredLength())
 
+    def testMaxNFractionAllPassNoNs(self):
+        """
+        Test filtering by maximum fraction of Ns. If there are no Ns in the
+        sequences, all must pass the filtering.
+        """
+        read1 = Read('id1', 'ATTA')
+        read2 = Read('id2', 'ATTAAC')
+        initialReads = Reads([read1, read2])
+        initialReads.filter(maxNFraction=0.9)
+
+        reads = Reads(initialReads)
+        self.assertEqual([read1, read2], list(reads))
+
+    def testMaxNFractionOnePasses(self):
+        """
+        Test filtering by maximum fraction of Ns. If there are too many Ns in
+        one of the sequences, only one must pass the filtering.
+        """
+        read1 = Read('id1', 'ATTA')
+        read2 = Read('id2', 'ATTNNN')
+        initialReads = Reads([read1, read2])
+        initialReads.filter(maxNFraction=0.4)
+
+        reads = Reads(initialReads)
+        self.assertEqual([read1], list(reads))
+
+    def testMaxNFractionAllPassNs(self):
+        """
+        Test filtering by maximum fraction of Ns. If there are Ns in the
+        sequence, but below the threshold, all sequences must pass the
+        filtering.
+        """
+        read1 = Read('id1', 'ATTA')
+        read2 = Read('id2', 'ATTNNN')
+        initialReads = Reads([read1, read2])
+        initialReads.filter(maxNFraction=0.6)
+
+        reads = Reads(initialReads)
+        self.assertEqual([read1, read2], list(reads))
+
     def testNoVariableSitesConfirm(self):
         """
         If two Reads have no bases that are variable, nothing should be