Converting domain-specific regular-expressions to a list of all matching instances

Question

There seem to be several questions floating around Stackexchange regarding how to take a python regular expression list the matching instances. This problem is a bit different because 1) I'm need to start off with a domain-specific regex system, not just the standard python regex, and 2) the possible matches for any given location in string are limited to a specific 'alphabet' of characters.

The current implementation works, but I've got a nagging feeling there's a much more succinct way to do it.

import numpy as np

def shorthand_to_instances(shorthand,alphabet):
    """
    Converts the shorthand-notation for bioinformatics 'motifs'
    into the all the instances of sequences they would match.
    Brackets means "it could be any of these letters" in this location.
    Curly braces means "it's *not* one of these letters" in this location,
    it can be any letter in the given "alphabet" except these.
    Lowercase "x" means it could be any of the letters in the given
    "alphabet".

    >>> binf_motifs_to_sequences('T[AG]','AGCT')
    ['TA', 'TG']

    >>> binf_motifs_to_sequences('T{A}T','AGCT')
    ['TGT', 'TCT', 'TTT']

    >>> binf_motifs_to_sequences('G{AGC}AGTC','AGCT')
    ['GTAGTC']

    >>> binf_motifs_to_sequences('xGTC','AGCT')
    ['AGTC', 'GGTC', 'CGTC', 'TGTC']

    """    
    instances=[]
    shc = shorthand
    while len(shc)>0:
        char = shc[0]
        assert char in alphabet+'['+']'+'{'+'}'+'x', "Failed assert: '{0}'".format(char)
        if char in alphabet:
            # it's in the alphabet
            if len(instances)==0:
                instances.append(char)
            else:
                for i in list(np.arange(len(instances))):
                    instances[i] += char
            shc = shc[1:]
        elif char=='[':
            # it's a 'one of these'            
            chargroup,c,e = shc[1:].partition(']')
            if len(instances)==0:
                for char in chargroup:
                    instances.append(char)
            else:
                newinstances = []
                for i in list(np.arange(len(instances))):
                    stem = instances[i]
                    for char in chargroup:
                        newinstances.append(stem+char)              
                instances = newinstances
            shc = shc[len(chargroup)+2:]
        elif char=='{':
            # it's a 'not this / not these'
            chargroup0,c,e = shc[1:].partition('}')
            chargroup = alphabet.translate(None,chargroup0)
            if len(instances)==0:
                for char in chargroup:
                    instances.append(char)
            else:
                newinstances = []
                for i in list(np.arange(len(instances))):
                    stem = instances[i]
                    for char in chargroup:
                        newinstances.append(stem+char)
                instances = newinstances
            shc = shc[len(chargroup0)+2:]
        elif char=='x':
            # it's an 'any char'
            chargroup = alphabet
            if len(instances)==0:
                for char in chargroup:
                    instances.append(char)
            else:
                newinstances = []
                for i in list(np.arange(len(instances))):
                    stem = instances[i]
                    for char in chargroup:
                        newinstances.append(stem+char)
                instances = newinstances
            shc = shc[1:]
        else:
            break

    return instances

if __name__ == "__main__":
    import doctest
    doctest.testmod()

Answer 1

1. Review

1. This is clear, well-documented code.

2. The doctests all call binf_motifs_to_sequences but your function is named shorthand_to_instances so all the doctests fail.

3. There's no need to have this kind of infrastructure for running doctests:

   if __name__ == "__main__":
       import doctest
       doctest.testmod()
   

   since you can easily run doctests from the command line using Python's -m option:

   python -m doctest mymodule.py
   

4. Assertions should be used for programming errors (conditions that you believe can't happen), but this assertion:

   assert char in alphabet+'['+']'+'{'+'}'+'x'
   

   is really a runtime error (it might happen if bad data was passed to the function), so it would be better to raise an exception.

5. A given shorthand can have exponentially many matching sequences. For example, xxxx already has 256 matching sequences over a four-character alphabet. For this reason, I would avoid returning a list of matching sequences, as this requires building a list and storing all the sequences in memory. It would be better to generate the sequences one at a time.

2. Simplifying

The way to make this kind of code simpler is to use functional decomposition to divide it into pieces that are simpler to understand, implement, and test.

A sensible division into two pieces would be to first, parse the shorthand notation, producing some kind of data structure; and second, take the data structure and generate the matching sequences. Since the alphabet is likely to be small, the natural data structure to use is just a sequence of strings giving the characters that match at each position. So given the shorthand xT[AG]{AT} over the alphabet ACGT, the intermediate data would be the sequence ['ACGT', 'T', 'AG', 'CG'].

(If alphabets were expected to be much larger, then this might be somewhat wasteful in space, and so there might be a need to have a more sophisticated data structure here.)

I've called the first step "parsing" but really it's so simple that it's just lexical analysis and in Python this can be done simply by writing a regular expression that matches a token, and calling re.finditer to loop over the tokens:

import re

class ParseError(Exception): pass

def shorthand_parse(shorthand, alphabet):
    """Parse bioinformatics shorthand notation and generate strings giving
    the characters in alphabet that are allowed to match at each
    position in the sequence.

    Square brackets means "it could be any of these letters" in this
    location. Curly braces means "it's *not* one of these letters" in
    this location (it can be any letter in the given alphabet except
    these). Lowercase "x" means it could be any of the letters in the
    given alphabet.

    >>> list(shorthand_parse('xT[AG]{AT}', 'ACGT'))
    ['ACGT', 'T', 'AG', 'CG']

    """
    # Use set() to remove duplicates, and sorted() to ensure that the
    # output is in a canonical order (for predictability and testing).
    alphabet = set(alphabet)
    alphabet_str = ''.join(sorted(alphabet))    
    shorthand_re = r'''(x)             # Any character in alphabet
                      |([{0}])         # Single character
                      |\[([{0}]*)\]    # One of these characters
                      |{{([{0}]*)}}    # Not one of these characters
                      |(.)             # Anything else is an error
                    '''.format(alphabet_str)
    for m in re.finditer(shorthand_re, shorthand, re.VERBOSE):
        x, single, oneof, noneof, error = m.groups()
        if x:
            yield alphabet_str
        elif single:
            yield single
        elif oneof:
            yield ''.join(sorted(set(oneof)))
        elif noneof:
            yield ''.join(sorted(alphabet - set(noneof)))
        else:
            raise ParseError("unexpected character '{}'".format(error))

Now that we have the strings of characters that are valid at each position, we want the Cartesian product of those strings. This can be computed by itertools.product, and in fact the output of shorthand_parse is already in the right format for input to that function:

import itertools

def shorthand_sequences(shorthand, alphabet):
    """Generate all sequences matching the given bioinformatics shorthand
    notation over the given alphabet.

    >>> list(shorthand_sequences('T[AG]', 'AGCT'))
    ['TA', 'TG']
    >>> list(shorthand_sequences('T{A}T', 'AGCT'))
    ['TCT', 'TGT', 'TTT']
    >>> list(shorthand_sequences('G{AGC}AGTC', 'AGCT'))
    ['GTAGTC']
    >>> list(shorthand_sequences('xGTC', 'AGCT'))
    ['AGTC', 'CGTC', 'GGTC', 'TGTC']
    >>> list(shorthand_sequences('[AG]{AT}[CG]', 'AGCT'))
    ['ACC', 'ACG', 'AGC', 'AGG', 'GCC', 'GCG', 'GGC', 'GGG']

    """    
    for seq in itertools.product(*shorthand_parse(shorthand, alphabet)):
        yield ''.join(seq)

Answer 2

I think an alternative approach that would be helpful is to separate out the parsing from generating the matching sequences. Under one method of parsing, the generation code is very succinct.

As an example take the call:

binf_motifs_to_sequences('xGTC[TA]','AGCT')

A parsed version of the shorthand might produce the following list:

[
 'AGCT', # 'x', any character from the alphabet
 'GTC',
 'TA'    # the bracketed characters
]

Taking the cartesian product of all the elements of this list produces the matching sequences. In python2.7 this looks like:

from itertools import product
product(*parsing_output)

---

---

I took a stab at bringing these concepts together in the code below. The parsing function is a simple left to right parse. It does not support any nested brackets/curlies. Admittedly, it's not that pretty. Of the course with separating out the parsing as a separate concern, it's easier to plugin something cleaner. The getMatchingSequences() function returns a generator instead of a list, which is against your specification. I figured that the conversion to a list is simple, and that it allows for space efficiency if using complex shorthands.

import itertools

class Motif:
    def __init__(self, shorthand, alphabet):
        self.shorthand = shorthand
        self.alphabet = alphabet
        self.special_chars = ['x', '[', ']', '{', '}']

        if not self.hasValidCharacters():
            raise ValueError('Invalid character(s) in shorthand %(shorthand)s' %
                               {'shorthand': self.shorthand})

    def hasValidCharacters(self):
        entire_alphabet = set(self.special_chars).union(set(self.alphabet))
        return set(self.shorthand).issubset(entire_alphabet)

    def parse(self):
        parsed = []
        curr_sequence = []
        i = 0
        while i < len(self.shorthand):
            char = self.shorthand[i]

            if char not in self.alphabet:
                # then is a special character
                if len(curr_sequence) > 0:
                    parsed.append(curr_sequence)
                    curr_sequence = []

                if char == 'x':
                    parsed.append(self.alphabet)
                elif char == '[':
                    i += 1
                    try:
                        right_brack_index = self.shorthand[i:].index(']') + i
                    except ValueError as e:
                        raise ValueError('missing right bracket "]" in shorthand: %(shorthand)s' %
                                {'shorthand': self.shorthand})

                    bracket_group = self.shorthand[i:right_brack_index]
                    i = right_brack_index
                    parsed.append(bracket_group)
                elif char == '{':
                    i += 1
                    try:
                        right_curly_index = self.shorthand[i:].index('}') + i
                    except ValueError as e:
                        raise ValueError('missing right curly "}" in shorthand: %(shorthand)s' %
                                {'shorthand': self.shorthand})

                    curly_group = self.shorthand[i:right_curly_index]
                    i = right_curly_index
                    alphabet_diff = set(self.alphabet) - set(curly_group)
                    parsed.append(list(alphabet_diff))
                elif char == '}' or ']':
                    raise ValueError('"%(char)s" appears without its lefthand equivalent' %
                                    {'char': char})
            else:
                if i < len(self.shorthand):
                    curr_sequence += self.shorthand[i]
            i += 1

        return parsed

    def getMatchingSequences(self):
        for sequence in itertools.product(*self.parse()):
            yield ''.join(sequence)

motif = Motif('xGTC[TA]', 'AGCT')
print list(motif.getMatchingSequences())