Repository: acorg/seqgen

Issue:
Add a 'linear' directive to allow a set of mutant sequences to be successively generated from the previous one

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index 4234921..29b9dee 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,4 @@ __pycache__
 build
 dist
 seqgen.egg-info
+.pytest_cache
diff --git a/README.md b/README.md
index 2ae492c..af241af 100644
--- a/README.md
+++ b/README.md
@@ -50,7 +50,7 @@ optional arguments:
                         100)
 ```
 
-## Sequence spefication
+## Sequence specification
 
 Your JSON specifies what sequences you want created.
 
@@ -148,6 +148,11 @@ specification keys you can put in a object in the JSON is as follows:
   an internal name.
 * `random aa`: The sequence should be made of random amino acids.
 * `random nt`: The sequence should be made of random nucleotides.
+* `ratchet`: If `true` and a count and mutation rate are given, sequences
+  after the first will be mutants of the previous sequence. This can be
+  used to build a series of sequences that are successive mutants of one
+  another. For the name derivation, see
+  [Muller's Ratchet](https://en.wikipedia.org/wiki/Muller's_ratchet).
 * `sections`: Gives a list of sequence sections used to build up another
   sequence (see example above). Each section is a JSON object that may
   contain the keys present in a regular sequence object, excluding `count`,
diff --git a/seqgen/__init__.py b/seqgen/__init__.py
index 58c2ae3..285814d 100644
--- a/seqgen/__init__.py
+++ b/seqgen/__init__.py
@@ -2,6 +2,6 @@ from seqgen.sequences import Sequences
 
 # Note that the version string must have the following format, otherwise it
 # will not be found by the version() function in ../setup.py
-__version__ = '1.0.10'
+__version__ = '1.0.11'
 
 __all__ = ['Sequences']
diff --git a/seqgen/sequences.py b/seqgen/sequences.py
index 7e31515..3af0d90 100644
--- a/seqgen/sequences.py
+++ b/seqgen/sequences.py
@@ -23,6 +23,7 @@ class Sequences(object):
     DEFAULT_LENGTH = 100
     DEFAULT_ID_PREFIX = 'seq-id-'
     LEGAL_SPEC_KEYS = {
+        'alphabet',
         'count',
         'description',
         'id',
@@ -33,12 +34,14 @@ class Sequences(object):
         'name',
         'random aa',
         'random nt',
+        'ratchet',
         'sections',
         'sequence',
         'sequence file',
         'skip',
     }
     LEGAL_SPEC_SECTION_KEYS = {
+        'alphabet',
         'from name',
         'length',
         'mutation rate',
@@ -82,8 +85,30 @@ class Sequences(object):
         self._vars = vars_
         self._sequenceSpecs = list(map(self._expandSpec, sequenceSpecs))
         self._checkKeys()
+        self._checkValid()
         self._names = {}
 
+    def _checkValid(self):
+        """
+        Check that all specification dicts contain sensible values.
+
+        @param sequenceSpec: A C{dict} with information about the sequences
+            to be produced.
+        @raise ValueError: If any problem is found.
+        """
+        for specCount, spec in enumerate(self._sequenceSpecs, start=1):
+            if spec.get('ratchet'):
+                nSequences = spec.get('count', 1)
+                if nSequences == 1:
+                    raise ValueError(
+                        'Sequence specification %d is specified as ratchet '
+                        'but its count is only 1.' % specCount)
+
+                if 'mutation rate' not in spec:
+                    raise ValueError(
+                        'Sequence specification %d is specified as ratchet '
+                        'but does not give a mutation rate.' % specCount)
+
     def _checkKeys(self):
         """
         Check that all specification dicts only contain legal keys.
@@ -147,11 +172,15 @@ class Sequences(object):
             new[k] = value
         return new
 
-    def _specToRead(self, spec):
+    def _specToRead(self, spec, previousRead=None):
         """
         Get a sequence from a specification.
 
         @param spec: A C{dict} with keys/values specifying a sequence.
+        @param previousRead: If not C{None}, a {dark.Read} instance containing
+            the last read this method returned. This is only used when
+            'ratchet' is given for a specification, in which case we generate
+            a mutant based on the previous read.
         @raise ValueError: If the section spec refers to a non-existent other
             sequence, or to part of another sequence but the requested part
             exceeds the bounds of the other sequence. Or if the C{spec} does
@@ -159,10 +188,14 @@ class Sequences(object):
             to.
         @return: A C{dark.Read} instance.
         """
-        nt = True
+        alphabet = self.NT
         length = spec.get('length', self._defaultLength)
 
-        if 'from name' in spec:
+        if spec.get('ratchet') and previousRead:
+            read = Read(None, previousRead.sequence)
+            alphabet = previousRead.alphabet
+
+        elif 'from name' in spec:
             name = spec['from name']
             try:
                 namedRead = self._names[name]
@@ -176,6 +209,7 @@ class Sequences(object):
                 # named read.
                 length = spec.get('length', len(namedRead))
                 sequence = namedRead.sequence[index:index + length]
+                alphabet = namedRead.alphabet
 
                 if len(sequence) != length:
                     raise ValueError(
@@ -201,37 +235,39 @@ class Sequences(object):
                 raise ValueError("Sequence file '%s' could not be read." %
                                  spec['sequence file'])
 
+        elif spec.get('alphabet'):
+            alphabet = spec['alphabet']
+            read = Read(None, ''.join(choice(alphabet) for _ in range(length)))
+
         elif spec.get('random aa'):
-            read = Read(None, ''.join(choice(self.AA) for _ in range(length)))
-            nt = False
+            alphabet = self.AA
+            read = Read(None, ''.join(choice(alphabet) for _ in range(length)))
 
         else:
-            # Assume random nucleotides are wanted.
-            read = Read(None, ''.join(choice(self.NT) for _ in range(length)))
+            read = Read(None, ''.join(choice(alphabet) for _ in range(length)))
 
         try:
             rate = spec['mutation rate']
         except KeyError:
             pass
         else:
-            read.sequence = self._mutate(read.sequence, rate, nt)
+            read.sequence = self._mutate(read.sequence, rate, alphabet)
+
+        read.alphabet = alphabet
 
         return read
 
-    def _mutate(self, sequence, rate, nt):
+    def _mutate(self, sequence, rate, alphabet):
         """
         Mutate a sequence at a certain rate.
 
         @param sequence: A C{str} nucleotide or amino acid sequence.
         @param rate: A C{float} mutation rate.
-        @param nt: If C{True} the sequence is nucleotides, else amino acids.
+        @param alphabet: A C{list} of alphabet letters.
         @return: The mutatated C{str} sequence.
         """
         result = []
-        if nt:
-            possibles = set(self.NT)
-        else:
-            possibles = set(self.AA)
+        possibles = set(alphabet)
         for current in sequence:
             if uniform(0.0, 1.0) < rate:
                 result.append(choice(list(possibles - {current})))
@@ -247,17 +283,24 @@ class Sequences(object):
         @param sequenceSpec: A C{dict} with information about the sequences
             to be produced.
         """
+        alphabet = None
+        previousRead = None
         nSequences = spec.get('count', 1)
         for count in range(nSequences):
             id_ = None
             if 'sections' in spec:
                 sequence = ''
                 for section in spec['sections']:
-                    sequence += self._specToRead(section).sequence
+                    read = self._specToRead(
+                        section, previousRead)
+                    sequence += read.sequence
+                    if alphabet is None:
+                        alphabet = read.alphabet
             else:
-                read = self._specToRead(spec)
+                read = self._specToRead(spec, previousRead)
                 sequence = read.sequence
                 id_ = read.id
+                alphabet = read.alphabet
 
             if id_ is None:
                 try:
@@ -284,6 +327,7 @@ class Sequences(object):
                 pass
 
             read = Read(id_, sequence)
+            read.alphabet = alphabet
 
             # Keep a reference to this result if it is named.
             if count == 0 and 'name' in spec:
@@ -296,6 +340,7 @@ class Sequences(object):
 
             if not spec.get('skip'):
                 yield read
+                previousRead = read
 
     def __iter__(self):
         for sequenceSpec in self._sequenceSpecs:


Test patch:
diff --git a/test/testSequences.py b/test/testSequences.py
index 5e9e4ef..b4b839d 100644
--- a/test/testSequences.py
+++ b/test/testSequences.py
@@ -33,7 +33,7 @@ class TestSequences(TestCase):
         If the specification is not valid JSON, a ValueError must be raised.
         """
         if PY3:
-            error = '^Expecting value: line 1 column 1 \(char 0\)$'
+            error = '^Expecting value: line 1 column 1 \\(char 0\\)$'
         else:
             error = '^No JSON object could be decoded$'
         assertRaisesRegex(self, ValueError, error, Sequences,
@@ -45,7 +45,7 @@ class TestSequences(TestCase):
         If the JSON specification has no 'sequences' key, a ValueError
         must be raised.
         """
-        error = "^The specification JSON must have a 'sequences' key\.$"
+        error = "^The specification JSON must have a 'sequences' key\\.$"
         assertRaisesRegex(self, ValueError, error, Sequences,
                           spec='filename')
 
@@ -57,7 +57,7 @@ class TestSequences(TestCase):
         must be raised.
         """
         s = Sequences(spec='filename')
-        error = "^Name 'a' is duplicated in the JSON specification\.$"
+        error = "^Name 'a' is duplicated in the JSON specification\\.$"
         assertRaisesRegex(self, ValueError, error, list, s)
 
     def testNoSequences(self):
@@ -137,18 +137,53 @@ class TestSequences(TestCase):
                 }
             ]
         }'''))
-        error = ("^Sequence with id 'the-id' has a count of 6\. If you want "
-                 "to specify one sequence with an id, the count must be 1\. "
+        error = ("^Sequence with id 'the-id' has a count of 6\\. If you want "
+                 "to specify one sequence with an id, the count must be 1\\. "
                  "To specify multiple sequences with an id prefix, use "
-                 "'id prefix'\.$")
+                 "'id prefix'\\.$")
         assertRaisesRegex(self, ValueError, error, list, s)
 
+    def testRatchetWithNoCount(self):
+        """
+        If ratchet is specified for seqeunce spec its count must be greater
+        than one.
+        """
+        spec = StringIO('''{
+            "sequences": [
+                {
+                    "id": "xxx",
+                    "ratchet": true
+                }
+            ]
+        }''')
+        error = ("^Sequence specification 1 is specified as ratchet but its "
+                 "count is only 1\\.$")
+        assertRaisesRegex(self, ValueError, error, Sequences, spec)
+
+    def testRatchetWithNoMutationRate(self):
+        """
+        If ratchet is specified for seqeunce spec it must have a mutation
+        rate.
+        """
+        spec = StringIO('''{
+            "sequences": [
+                {
+                    "count": 4,
+                    "id": "xxx",
+                    "ratchet": true
+                }
+            ]
+        }''')
+        error = ("^Sequence specification 1 is specified as ratchet but does "
+                 "not give a mutation rate\\.$")
+        assertRaisesRegex(self, ValueError, error, Sequences, spec)
+
     def testUnknownSpecificationKey(self):
         """
         If an unknown key is given in a sequence specification, a ValueError
         must be raised.
         """
-        error = "^Sequence specification 1 contains an unknown key: dog\.$"
+        error = "^Sequence specification 1 contains an unknown key: dog\\.$"
         assertRaisesRegex(self, ValueError, error, Sequences, StringIO('''{
             "sequences": [
                 {
@@ -164,7 +199,7 @@ class TestSequences(TestCase):
         """
         for key in 'id', 'id prefix', 'name':
             error = ("^Section 1 of sequence specification 1 contains an "
-                     "unknown key: %s\.$" % key)
+                     "unknown key: %s\\.$" % key)
             assertRaisesRegex(self, ValueError, error, Sequences, StringIO('''{
                 "sequences": [
                     {
@@ -177,6 +212,34 @@ class TestSequences(TestCase):
                 ]
             }''' % key))
 
+    def testOneLetterAlphabet(self):
+        """
+        It must be possible to specify an alphabet with just one symbol.
+        """
+        s = Sequences(StringIO('''{
+            "sequences": [
+                {
+                    "alphabet": "0"
+                }
+            ]
+        }'''), defaultLength=500)
+        (read,) = list(s)
+        self.assertEqual('0' * 500, read.sequence)
+
+    def testTwoLetterAlphabet(self):
+        """
+        It must be possible to specify an alphabet with two symbols.
+        """
+        s = Sequences(StringIO('''{
+            "sequences": [
+                {
+                    "alphabet": "01"
+                }
+            ]
+        }'''), defaultLength=500)
+        (read,) = list(s)
+        self.assertTrue(x in '01' for x in read.sequence)
+
     def testOneSequenceIdOnlyDefaultLength(self):
         """
         If only one sequence is specified, and only by id, one sequence
@@ -526,7 +589,7 @@ class TestSequences(TestCase):
             ]
         }'''))
         error = ("^Sequence section refers to name 'xxx' of "
-                 "non-existent other sequence\.$")
+                 "non-existent other sequence\\.$")
         assertRaisesRegex(self, ValueError, error, list, s)
 
     def testSectionWithNameReferenceTooShort(self):
@@ -552,7 +615,7 @@ class TestSequences(TestCase):
         }'''))
         error = ("^Sequence specification refers to sequence name 'xxx', "
                  "starting at index 1 with length 10, but 'xxx' is not long "
-                 "enough to support that\.$")
+                 "enough to support that\\.$")
         assertRaisesRegex(self, ValueError, error, list, s)
 
     def testNamedRecombinant(self):
@@ -636,3 +699,42 @@ class TestSequences(TestCase):
         diffs = sum((a != b) for (a, b) in zip(sequence, read.sequence))
         self.assertEqual(len(sequence), len(read.sequence))
         self.assertEqual(diffs, len(read.sequence))
+
+    def testRatchet(self):
+        """
+        The ratchet specification must result in the expected result.
+        """
+        # Note that this is a very simple test, using a 1.0 mutation rate
+        # and a fixed alphabet.
+        length = 50
+        s = Sequences(StringIO('''{
+            "sequences": [
+                {
+                    "name": "orig",
+                    "alphabet": "01",
+                    "length": %s
+                },
+                {
+                    "count": 2,
+                    "from name": "orig",
+                    "mutation rate": 1.0,
+                    "ratchet": true
+                }
+            ]
+        }''' % length))
+        (orig, mutant1, mutant2) = list(s)
+        # The distance from the original to the first mutant must be 100 (i.e.,
+        # all bases).
+        diffCount = sum(a != b for (a, b) in
+                        zip(orig.sequence, mutant1.sequence))
+        self.assertEqual(length, diffCount)
+
+        # The distance from the first mutant to the second must be 100 (i.e.,
+        # all bases).
+        diffCount = sum(a != b for (a, b) in
+                        zip(mutant1.sequence, mutant2.sequence))
+        self.assertEqual(length, diffCount)
+
+        # The sequences of the original and the second mutant must be
+        # identical.
+        self.assertEqual(orig.sequence, mutant2.sequence)