13
\$\begingroup\$

The goal with this function is to find one DNA sequence within another sequence, with a specified amount of mismatch tolerance. For example:

  • dnasearch('acc','ccc',0) shouldn't find a match, while
  • dnasearch('acc','ccc',1) should find one.

EDIT: Function should also be able to find pat sequence in a sequence larger than the pat sequence; furthermore it should be able to find multiple matches within that sequence if they exist. For example, dnasearch('acc','acaaatc',1) would find two matches: 'aca' and 'act'. Thanks to @weronika for pointing this out.

Here's my working draft:

def dnasearch(pat,seq,mismatch_allowed=0):
patset1=set([pat])
patset2=set([pat])
for number in range(mismatch_allowed):
    for version in patset1:
        for a in range(len(version)):
            for letter in 'actg':
                newpat=version[:a]+letter+version[a+1:]
                patset2.add(newpat)
    patset1|=patset2
for n in patset1:
    if seq.find(n) != -1: print 'At position', seq.find(n),\
       'found',n

I was wondering how to write this function:

  1. With fewer for-loops.
  2. Without the second patset. Originally I tried just adding new patterns to patset1, but I got the error: RuntimeError: Set changed size during iteration.
  3. Simpler in any other way.
\$\endgroup\$
4
  • 1
    \$\begingroup\$ Do you mean en.wikipedia.org/wiki/Levenshtein_distance (insertion, deletion, or substitution of a single character)? \$\endgroup\$
    – cat_baxter
    Jun 13, 2012 at 8:42
  • \$\begingroup\$ You need to add a few more examples to clarify the problem, I think. It's unclear from the current examples what should happen when the two input sequences are different lengths. \$\endgroup\$
    – weronika
    Jun 13, 2012 at 19:01
  • \$\begingroup\$ Look like InterviewStreet puzzle "Save Humanity" :) \$\endgroup\$
    – cat_baxter
    Jun 14, 2012 at 10:20
  • \$\begingroup\$ 2012: where the words "simple" and "DNS sequence finder" appear in the same sentence. \$\endgroup\$
    – msanford
    Jun 14, 2012 at 19:34

1 Answer 1

4
\$\begingroup\$

This is a solved problem in bioinformatics and there are specific algorithms depending on your specific use-case:

  • Do you allow gaps?
  • How many mismatches do gaps incur?
  • … etc.

In the general case, this is a global pairwise alignment which is computed by the Needleman–Wunsch algorithm. If you are just interested in the score of the alignment, rather than the complete alignment, this can be simplified by not performing the backtracking.

Now, if you have a threshold and are only interested in finding matches below said threshold, then a more efficient variant employs the Ukkonen trick – unfortunately, references for this are hard to find online, and the print literature is quite expensive. But what it does is essentially break the recurrence in the above-mentioned algorithm once the score exceeds the threshold chosen before (using the insight that the error score can only increase, never decrease in a given column).

But all this is unnecessary if you don’t allow gaps. In that case, the algorithm is as straightforward as walking over both strings at once and looking whether the current positions match.

And this can be expressed in a single line:

def distance(a, b):
    return sum(map(lambda (x, y): 0 if x == y else 1, zip(a, b)))

and:

def distance_less_than(k, a, b):
    return distance(a, b) < k
\$\endgroup\$
5
  • \$\begingroup\$ Thanks for the answer Konrad--I love the distance function in your answer. I did forget an important detail in my question though, pointed out by weronika, which is that the function should also be able to find the pattern sequence in a sequence larger than itself. And yes, I'm not worrying about gaps for the time being. \$\endgroup\$
    – Sean
    Jun 13, 2012 at 19:24
  • \$\begingroup\$ The distance function can be written more simply def distance(a, b): return sum(x == y for x, y in zip(a, b))) \$\endgroup\$
    – Caridorc
    Apr 7, 2015 at 12:30
  • \$\begingroup\$ @Caridorc That’s indeed correct. However, I prefer not muddling boolean and numeric types. In fact, I wish Python wouldn’t allow your solution, because it relies on weak typing. \$\endgroup\$ Apr 7, 2015 at 14:08
  • \$\begingroup\$ That is impossibile because Python is strongly typed, bool is a subset of int. \$\endgroup\$
    – Caridorc
    Apr 7, 2015 at 15:14
  • 1
    \$\begingroup\$ @Caridorc “strongly typed” is a relative, not an absolute statement. Python is indeed relatively strongly typed (compared to other languages), which is why I have zero tolerance for this “a bool is an integer” nonsense. Logically, a boolean is a quite distinct concept from a nominally unbounded integer number, and it rarely makes sense to treat one as a subtype of the other (in fact, it leads to a ton of nonsense; case in point: if we treat booleans as numbers then they are obviously a ring modulo 2 and True+True should either be False or True, not 2, because of modular arithmetic rules) \$\endgroup\$ Apr 7, 2015 at 15:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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