12
\$\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\$
  • 1
    \$\begingroup\$ Do you mean en.wikipedia.org/wiki/Levenshtein_distance (insertion, deletion, or substitution of a single character)? \$\endgroup\$ – cat_baxter Jun 13 '12 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 '12 at 19:01
  • \$\begingroup\$ Look like InterviewStreet puzzle "Save Humanity" :) \$\endgroup\$ – cat_baxter Jun 14 '12 at 10:20
  • \$\begingroup\$ 2012: where the words "simple" and "DNS sequence finder" appear in the same sentence. \$\endgroup\$ – msanford Jun 14 '12 at 19:34
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\$
  • \$\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 '12 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 '15 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\$ – Konrad Rudolph Apr 7 '15 at 14:08
  • \$\begingroup\$ That is impossibile because Python is strongly typed, bool is a subset of int. \$\endgroup\$ – Caridorc Apr 7 '15 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\$ – Konrad Rudolph Apr 7 '15 at 15:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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