Finding Subsequences

To demonstrate how to find the subsequence, I’ll start by defining the following sequence and subsequence in the REPL:

>>> seq = 'GATATATGCATATACTT'
>>> subseq = 'ATAT'

I can use in to determine if one sequence is a subset of another. This also works for membership in lists, sets, or keys of a dictionary:

>>> subseq in seq
True

That’s good information, but it doesn’t tell me where the string can be found. Luckily there’s the str.find() function, which says subseq can be found starting at index 1 (which is the second character):

>>> seq.find(subseq)
1

I know from the Rosalind description that the answer should be 2, 4, and 10. I just found 2, so how can I find the next? I can’t just call the same function again because I’ll get the same answer. I need to look further into the sequence. Maybe help(str.find) could be of some use?

>>> help(str.find)
find(...)
    S.find(sub[, start[, end]]) -> int

    Return the lowest index in S where substring sub is found,
    such that sub is contained within S[start:end].  Optional
    arguments start and end are interpreted as in slice notation.

    Return -1 on failure.

It appears I can specify a start position. I’ll use 1 greater than the position where the first subsequence was found, which was 1, so starting at 2:

>>> seq.find(subseq, 2)
3

Great. That’s the next answer—well, 4 is the next answer, but you know what I mean. I’ll try that again, this time starting at 4:

>>> seq.find(subseq, 4)
9

That was the last value I expected. What happens if I try using a start of 10? As the documentation shows, this will return -1 to indicate the subsequence cannot be found:

>>> seq.find(subseq, 10)
-1

Can you think of a way to iterate through the sequence, remembering the last position where the subsequence was found until it cannot be found?

Another option would be to use str.index(), but only if the subsequence is present:

>>> if subseq in seq:
...     seq.index(subseq)
...
1

To find the next occurrence, you could slice the sequence using the last known position. You’ll have to add this position to the starting position, but you’re essentially doing the same operation of moving further into the sequence to find if the subsequence is present and where:

>>> if subseq in seq[2:]:
...     seq.index(subseq[2:])
...
1

If you read help(str.index), you’ll see that, like str.find(), the function takes a second optional start position of the index at which to start looking:

>>> if subseq in seq[2:]:
...     seq.index(subseq, 2)
...
3

A third approach would be to use k-mers. If the subsequence is present, then it is by definition a k-mer, where k is the length of the subsequence. Use your code from Chapter 7 to extract all the k-mers and their positions from the sequence, and note the positions of the k-mers that match the subsequence.

Finally, since I’m looking for a pattern of text, I could use a regular expression. In Chapter 5, I used the re.findall() function to find all the Gs and Cs in DNA. I can similarly use this method to find all the subsequences in the sequence:

>>> import re
>>> re.findall(subseq, seq)
['ATAT', 'ATAT']

That seems to have a couple of problems. One is that it only returned two of the subsequences when I know there are three. The other problem is that this provides no information about where the matches are found. Fear not, the re.finditer() function solves this second problem:

>>> list(re.finditer(subseq, seq))
[<re.Match object; span=(1, 5), match='ATAT'>,
 <re.Match object; span=(9, 13), match='ATAT'>]

Now it’s apparent that it finds the first and last subsequences. Why doesn’t it find the second instance? It turns out regular expressions don’t handle overlapping patterns very well, but some additions to the search pattern can fix this. I’ll leave it to you and some internet searching to see if you can figure out a solution.

I’ve presented four different options for how to solve this problem. See if you can write solutions using each approach. The point is to explore the corners of Python, storing away tasty bits and tricks that might prove decisive in some future program you write. It’s OK to spend hours or days working this out. Keep at it until you have solutions that pass both pytest and make test.

Solution 2: Using the str.index() Method

This next solution is a variation that slices the sequence using the last known position:

def main() -> None:
    args = get_args()
    seq, subseq = args.seq, args.subseq 
    found = []
    last = 0
    while subseq in seq[last:]: 
        last = seq.index(subseq, last) + 1 
        found.append(last) 

    print(' '.join(map(str, found))) 

Unpack the sequence and subsequence.

Ask if the subsequence appears in a slice of the sequence starting at the last found position. The while loop will execute as long as this condition is true.

Use str.index() to get the starting position of the subsequence. The last variable gets updated by adding 1 to the subsequence index to create the next starting position.

Append this position to the list of found positions.

Use map() to coerce all the found integer positions to strings, then join them on spaces to print.

Here again, I rely on tracking the last place a subsequence was found. I start at position 0, or the beginning of the string:

>>> last = 0
>>> if subseq in seq[last:]:
...     last = seq.index(subseq, last) + 1
...
>>> last
2

The while True loop in the first solution is a common way to start an infinite loop. Here, the while loop will only execute as long as the subsequence is found in the slice of the sequence, meaning I don’t have to manually decide when to break out of the loop:

>>> last = 0
>>> found = []
>>> while subseq in seq[last:]:
...     last = seq.index(subseq, last) + 1
...     found.append(last)
...
>>> found
[2, 4, 10]

The found positions, in this case, are a list of integer values. In the first solution, I used *found to splat the list and relied on print() to coerce the values to strings and join them on spaces. If instead I were to try to create a new string from found using str.join(), I would run into problems. The str.join() function joins many strings into a single string and so raises an exception when you give it nonstring values:

>>> ' '.join(found)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: sequence item 0: expected str instance, int found

I could use a list comprehension to turn each number n into a string using the str() function:

>>> ' '.join([str(n) for n in found])
'2 4 10'

This can also be written using a map():

>>> ' '.join(map(lambda n: str(n), found))
'2 4 10'

I can leave out the lambda entirely because the str() function expects a single argument, and map() will naturally supply each value from found as the argument to str(). This is my preferred way to turn a list of integers into a list of strings:

>>> ' '.join(map(str, found))
'2 4 10'

Solution 3: A Purely Functional Approach

This next solution combines many of the preceding ideas using a purely functional approach. To start, consider the while loops in the first two solutions used to append nonnegative values to the found list. Does that sound like something a list comprehension could do? The range of values to iterate includes all the positions n from 0 to the end of the sequence minus the length of the subsequence:

>>> r = range(len(seq) - len(subseq))
>>> [n for n in r]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]

A list comprehension can use these values with str.find() to search for the subsequence in the sequence starting at each position n. Starting at positions 0 and 1, the subsequence can be found at index 1. Starting at positions 2 and 3, the subsequence can be found at index 3. This continues until -1 indicates the subsequence is not present for those positions n:

>>> [seq.find(subseq, n) for n in r]
[1, 1, 3, 3, 9, 9, 9, 9, 9, 9, -1, -1, -1]

I only want the nonnegative values, so I use filter() to remove them:

>>> list(filter(lambda n: n >= 0, [seq.find(subseq, n) for n in r]))
[1, 1, 3, 3, 9, 9, 9, 9, 9, 9]

Which could also be written by reversing the comparison in the lambda:

>>> list(filter(lambda n: 0 <= n, [seq.find(subseq, n) for n in r]))
[1, 1, 3, 3, 9, 9, 9, 9, 9, 9]

I show you this because I’d like to use partial() with the operator.le() (less than or equal) function because I don’t like lambda expressions:

>>> from functools import partial
>>> import operator
>>> ok = partial(operator.le, 0)
>>> list(filter(ok, [seq.find(subseq, n) for n in r]))
[1, 1, 3, 3, 9, 9, 9, 9, 9, 9]

I’d like to change the list comprehension to a map():

>>> list(filter(ok, map(lambda n: seq.find(subseq, n), r)))
[1, 1, 3, 3, 9, 9, 9, 9, 9, 9]

but again I want to get rid of the lambda by using partial():

>>> find = partial(seq.find, subseq)
>>> list(filter(ok, map(find, r)))
[1, 1, 3, 3, 9, 9, 9, 9, 9, 9]

I can use set() to get a distinct list:

>>> set(filter(ok, map(find, r)))
{1, 3, 9}

These are almost the correct values, but they are the index positions, which are zero-based. I need the values one greater, so I can make a function to add 1 and apply this using a map():

>>> add1 = partial(operator.add, 1)
>>> list(map(add1, set(filter(ok, map(find, r)))))
[2, 4, 10]

In these limited examples, the results are properly sorted; however, one can never rely on the order of values from a set. I must use the sorted() function to ensure they are properly sorted numerically:

>>> sorted(map(add1, set(filter(ok, map(find, r)))))
[2, 4, 10]

Finally, I need to print these values, which still exist as a list of integers:

>>> print(sorted(map(add1, set(filter(ok, map(find, r))))))
[2, 4, 10]

That’s almost right. As in the first solution, I need to splat the results to get print() to see the individual elements:

>>> print(*sorted(map(add1, set(filter(ok, map(find, r))))))
2 4 10

That’s a lot of closing parentheses. This code is starting to look a little like Lisp. If you combine all these ideas, you wind up with the same answer as the imperative solution but now by combining only functions:

def main() -> None:
    args = get_args()
    seq, subseq = args.seq, args.subseq 
    r = list(range(len(seq) - len(subseq))) 
    ok = partial(operator.le, 0) 
    find = partial(seq.find, subseq) 
    add1 = partial(operator.add, 1) 
    print(*sorted(map(add1, set(filter(ok, map(find, r)))))) 

Unpack the sequence and subsequence.

Generate a range of numbers up to the length of the sequence less the length of the subsequence.

Create a partial ok() function that will return True if a given number is greater than or equal to 0.

Create a partial find() function that will look for the subsequence in the sequence when provided with a start parameter.

Create a partial add1() function that will return one greater than the argument.

Apply all the numbers from the range to the find() function, filter out negative values, make the result unique by using the set() function, add one to the values, and sort the numbers before printing.

This solution uses only pure functions and would be fairly easy to understand for a person with a background in the Haskell programming language. If it seems like a jumble of confusion to you, I’d encourage you to spend some time working in the REPL with each piece until you understand how all these functions fit together perfectly.

Solution 4: Using K-mers

I mentioned that you might try finding the answer using k-mers, which I showed in Chapter 7. If the subsequence exists in the sequence, then it must be a k-mer, where k equals the length of the subsequence:

>>> seq = 'GATATATGCATATACTT'
>>> subseq = 'ATAT'
>>> k = len(subseq)
>>> k
4

Here are all the 4-mers in the sequence:

>>> kmers = [seq[i:i + k] for i in range(len(seq) - k + 1)]
>>> kmers
['GATA', 'ATAT', 'TATA', 'ATAT', 'TATG', 'ATGC', 'TGCA', 'GCAT', 'CATA',
 'ATAT', 'TATA', 'ATAC', 'TACT', 'ACTT']

Here are the 4-mers that are the same as the subsequence I’m looking for:

>>> list(filter(lambda s: s == subseq, kmers))
['ATAT', 'ATAT', 'ATAT']

I need to know the positions as well as the k-mers. The enumerate() function will return both the index and value of all the elements in a sequence. Here are the first four:

>>> kmers = list(enumerate([seq[i:i + k] for i in range(len(seq) - k + 1)]))
>>> kmers[:4]
[(0, 'GATA'), (1, 'ATAT'), (2, 'TATA'), (3, 'ATAT')]

I can use this with filter(), but now the lambda is receiving a tuple of the index and value so I will need to look at the second field (which is at index 1):

>>> list(filter(lambda t: t[1] == subseq, kmers))
[(1, 'ATAT'), (3, 'ATAT'), (9, 'ATAT')]

I really only care about getting the index for the matching k-mers. I could rewrite this using a map() with an if expression to return the index position when it’s a match, and None otherwise:

>>> list(map(lambda t: t[0] if t[1] == subseq else None, kmers))
[None, 1, None, 3, None, None, None, None, None, 9, None, None, None, None]

I’m frustrated that the standard map() function can only pass a single value to the lambda. What I need is some way to splat the tuple, like *t, to turn it into two values. Luckily, I’ve studied the itertools module documentation and found the starmap() function, so named because it will add a star to the lambda argument to splat it. This allows me to unpack a tuple value like (0, 'GATA') into the variables i with the index value of 0 and kmer with the value 'GATA'. With these, I can compare the kmer to the subsequence and also add 1 to the index (i):

>>> from itertools import starmap
>>> list(starmap(lambda i, kmer: i + 1 if kmer == subseq else None, kmers))
[None, 2, None, 4, None, None, None, None, None, 10, None, None, None, None]

This probably seems like an odd choice until I show you that filter(), if passed None for the lambda, will use the truthiness of each value, so that None values will be excluded. Because this line of code is getting rather long, I’ll write the function f() for map() on a separate line:

>>> f = lambda i, kmer: i + 1 if kmer == subseq else None
>>> list(filter(None, starmap(f, kmers)))
[2, 4, 10]

I can express a k-mer solution using imperative techniques:

def main() -> None:
    args = get_args()
    seq, subseq = args.seq, args.subseq
    k = len(subseq) 
    kmers = [seq[i:i + k] for i in range(len(seq) - k + 1)] 
    found = [i + 1 for i, kmer in enumerate(kmers) if kmer == subseq] 
    print(*found) 

When looking for k-mers, k is the length of the subsequence.

Use a list comprehension to generate all the k-mers from the sequence.

Iterate through the index and value of all the k-mers, where the k-mer is equal to the subsequence. Return one greater than the index position.

Print the found positions.

I can also express these ideas using purely functional techniques. Note that mypy insists on a type annotation for the found variable:

def main() -> None:
    args = get_args()
    seq, subseq = args.seq, args.subseq
    k = len(subseq)
    kmers = enumerate(seq[i:i + k] for i in range(len(seq) - k + 1)) 
    found: Iterator[int] = filter( 
        None, starmap(lambda i, kmer: i + 1 if kmer == subseq else None, kmers))
    print(*found) 

Generate an enumerated list of the k-mers.

Select the positions of those k-mers equal to the subsequence.

Print the results.

I find the imperative version easier to read, but would recommend you use whichever you find most intuitive. Whichever solution you prefer, the interesting point is that k-mers can prove extremely useful in many situations, such as for partial sequence comparison.

Solution 5: Finding Overlapping Patterns Using Regular Expressions

To this point, I’ve been writing fairly complex solutions to find a pattern of characters inside a string. This is precisely the domain of regular expressions, and so it’s a bit silly to write manual solutions. I showed earlier in this chapter that the re.finditer() function does not find overlapping matches and so returns just two hits when I know there are three:

>>> import re
>>> list(re.finditer(subseq, seq))
[<re.Match object; span=(1, 5), match='ATAT'>,
 <re.Match object; span=(9, 13), match='ATAT'>]

The problem turns out to be that the regex engine consumes strings as they match. That is, once the engine matches the first ATAT, it starts searching again at the end of the match. The solution is to wrap the search pattern in a look-ahead assertion using the syntax ?=(<pattern>) so that the engine won’t consume the matched string. Note that this is a positive look-ahead; there are also negative look-ahead assertions as well as both positive and negative look-behind assertions.

So if the subsequence is ATAT, then I want the pattern to be ?=(ATAT). The problem now is that the regex engine won’t save the match—I’ve just told it to look for this pattern but haven’t told it to do anything with the text that is found. I need to further wrap the assertion in parentheses to create a capture group:

>>> list(re.finditer('(?=(ATAT))', 'GATATATGCATATACTT'))
[<re.Match object; span=(1, 1), match=''>,
 <re.Match object; span=(3, 3), match=''>,
 <re.Match object; span=(9, 9), match=''>]

I can use a list comprehension over this iterator to call the match.start() function on each of the re.Match objects, adding 1 to correct the position:

>>> [match.start() + 1 for match in re.finditer(f'(?=({subseq}))', seq)]
[2, 4, 10]

Here is the final solution that I would suggest as the best way to solve this problem:

def main() -> None:
    args = get_args()
    seq, subseq = args.seq, args.subseq
    print(*[m.start() + 1 for m in re.finditer(f'(?=({subseq}))', seq)])