up vote 3 down vote favorite

For a bioinformatic problem, I found myself wanting to create a list of strings within Hamming distance "k" of a reference sequence. I wanted to do so quickly and pythonically. I've implemented it in pure python and in cython, with and without type declarations. The time performance has been identical. (I also compared the compiled python version with an interpreted version defined in ipython, which also performed identically.)

cfi has been set as shorthand to chain.from_iterable in order to reduce the number of dot operators being used as in the following module-level import and definition:

from itertools import chain
cfi = chain.from_iterable


@cython.returns(list)
def PermuteMotifOnce(cython.str motif, set alphabet={"A", "C", "G", "T"}):
    """
    Gets all strings within hamming distance 1 of motif and returns it as a
    list.
    """
    return list(set(cfi([[
        motif[:pos] + alpha + motif[pos + 1:] for
        alpha in alphabet] for
                         pos in range(len(motif))])))


def PyPermuteMotifOnce(motif, alphabet={"A", "C", "G", "T"}):
    """
    Gets all strings within hamming distance 1 of motif and returns it as a
    list.
    """
    return list(set(cfi([[
        motif[:pos] + alpha + motif[pos + 1:] for
        alpha in alphabet] for
                         pos in range(len(motif))])))


@cython.returns(list)
def PermuteMotifN(cython.str motif, cython.long n=-1):
    assert n > 0
    cdef set workingSet
    cdef cython.long i
    workingSet = {motif}
    for i in range(n):
        workingSet = set(cfi(map(PermuteMotifOnce, workingSet)))
    return list(workingSet)


def PyPermuteMotifN(motif, n=-1):
    assert n > 0
    workingSet = {motif}
    for i in range(n):
        workingSet = set(cfi(map(PermuteMotifOnce, workingSet)))
    return list(workingSet)

Results:

motif = "ACCTACTGAACT"
%timeit -n 5 PermuteMotifN(motif, 6)
5 loops, best of 3: 6.93s per loop
%timeit -n 5 PyPermuteMotifN(motif, 6)
5 loops, best of 3: 6.81s per loop
%timeit -n 5000 PyPermuteMotifN(motif, 2)
5000 loops, best of 3: 589 microseconds per loop
%timeit -n 5000 PermuteMotifN(motif, 2)
5000 loops, best of 3: 645 microseconds per loop

Is it just me, or does the pure Python seem even faster than the Cython? Is there a lot of time lost in extra type checking?

share|improve this question
    
For speed use pypy, for C-like speed use R-python – Caridorc May 5 '15 at 19:33
    
@Caridorc - I am not planning on using pypy, as many libraries I need are not compatible with it (or at least not compatible with it at a reasonable speed.) Cf. Link to related answer – NoSeatbelts May 5 '15 at 19:44
up vote 3 down vote accepted

Your strategy here is to keep a working set and repeatedly expand it by adding strings within Hamming distance 1 of elements in the set. This involves a lot of duplicated work, because the same candidate strings are going to be generated many times over (and then thrown away because they are already members of the working set).

I find that I get a substantial speedup just by generating each string exactly once.

from itertools import chain, combinations, product

def hamming_circle(s, n, alphabet):
    """Generate strings over alphabet whose Hamming distance from s is
    exactly n.

    >>> sorted(hamming_circle('abc', 0, 'abc'))
    ['abc']
    >>> sorted(hamming_circle('abc', 1, 'abc'))
    ['aac', 'aba', 'abb', 'acc', 'bbc', 'cbc']
    >>> sorted(hamming_circle('aaa', 2, 'ab'))
    ['abb', 'bab', 'bba']

    """
    for positions in combinations(range(len(s)), n):
        for replacements in product(range(len(alphabet) - 1), repeat=n):
            cousin = list(s)
            for p, r in zip(positions, replacements):
                if cousin[p] == alphabet[r]:
                    cousin[p] = alphabet[-1]
                else:
                    cousin[p] = alphabet[r]
            yield ''.join(cousin)

def hamming_ball(s, n, alphabet):
    """Generate strings over alphabet whose Hamming distance from s is
    less than or equal to n.

    >>> sorted(hamming_ball('abc', 0, 'abc'))
    ['abc']
    >>> sorted(hamming_ball('abc', 1, 'abc'))
    ['aac', 'aba', 'abb', 'abc', 'acc', 'bbc', 'cbc']
    >>> sorted(hamming_ball('aaa', 2, 'ab'))
    ['aaa', 'aab', 'aba', 'abb', 'baa', 'bab', 'bba']

    """
    return chain.from_iterable(hamming_circle(s, i, alphabet)
                               for i in range(n + 1))

In pure Python, this runs about 4 times as fast as the code in the original post.

share|improve this answer
    
Excellent - this has given me a better appreciation for the use of itertools as well. – NoSeatbelts May 5 '15 at 20:29
    
Can we do better? I find hamming_ball still kind of slow. – zyxue Sep 19 '16 at 22:30
    
@zyxue: There are exponentially many strings in the Hamming ball, so generating them all is always going to be slow, whatever method is used. Can you find a way to solve your problem without needing to generate them all? – Gareth Rees Sep 20 '16 at 14:58
    
@GarethRees, thank you for your reply. No, I have to generate them all. – zyxue Sep 20 '16 at 16:37
    
@zyxue: Well then, that might make a good question for Code Review. Explain the problem you are trying to solve, and present the code you've written so far, and then maybe someone will find a faster approach. – Gareth Rees Sep 24 '16 at 16:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.