Genome string clump finding problem

Question

I am trying to solve a bioinformatics problems from a Stepic course.

The problem posed: find clumps of the same pattern within a longer genome.

Motivation: Identifying 3 occurrences of the same pattern within a window of length 500 within a genome allows biologists to find the "replication origin" of a bacterial genome. This replication origin is required for the bacteria to clone itself, and it is therefore significant.

Using Python, I have written a solution that works effectively for small, test data sets, but which is not efficient for real-life genomes. Below I describe the method used.

I've used two core functions.

patternOccurances

This function finds all the starting locations of a short pattern in a longer string.

def patternOccurances2(text="GATATATGCATATACTT", pattern="ATAT"):

  occurances = []
  position = -1

  position = config.genome.find(pattern, position+1)
  occurances.append(position)

  while position != -1:
    position = text.find(pattern, position+1)
    occurances.append(position)
  return occurances[:-1]

frequentPatterns

This second function finds all patterns of length k in a string that occur at least t times; or, if t is omitted, only the most frequent patterns of length k.

def frequentPatterns(text="", k=2,numberOfOccurances=None):

  #Frequent Words Problem: Find the frequent k-mers in a string.
  #     Input: A string text, and an integer k, number of occurances t
  #     Output: All most frequent k-mers in Text.

  words = {}
  frequentPatterns = []

  #create a dictionary of all patterns of length k in text with their respective frequencies
  for i in range(len(text)-k + 1): #iterate over the valid length of the string
    a = text[i:i+k]
    if a in words:
      words[a] = words[a] + 1
    else:
      words[a] = 1

  if numberOfOccurances == None: #gives only the most frequent strings of length k

    largestValue=0
    for k, v in words.iteritems():
      if v > largestValue:
    frequentPatterns = []
    frequentPatterns.append(k) 
    largestValue = v
      if v == largestValue:
    frequentPatterns.append(k)

  else: #gives all strings of length k that occur numberOfOccurances times

    for k, v in words.iteritems():
      if v >= numberOfOccurances:
    frequentPatterns.append(k)

  return frequentPatterns

Main Code

The main code calls the above functions as follows:

#Clump Finding Problem: Find patterns forming clumps in a string.
#     Input: A long string Genome, and integers k, L, and t.
#            k is the length of the pattern we wish to match.
#            t is the minimum number of times we wish to find this pattern in a portion of length L in the Genome.
#     Output: All distinct k-mers forming (L, t)-clumps in Genome.


f = open("E-coli.txt",'r') #input from a large file, 4.2MB

genome =  f.read()

k = 9
L = 500
t = 3

#find all patterns of length k in genome that occur at least t times
frequentPatterns = ex213.frequentPatterns(text=genome,k=k,numberOfOccurances=t)

#for each pattern in frequentPatterns, find the locations of the patterns in the genome
#then find if any 3 of these patterns are within distance L of each other
for pattern in frequentPatterns:
  locations = ex322.patternOccurances("",pattern)
  for location in range(len(locations) - 2):
    if locations[location + 1] - locations[location] <= L:
      if locations[location + 2] - locations[location] <= L:
    print pattern #3 of the pattern are sufficiently close to each other
    break

f.close()

Main Code Performance

Measurements show that the frequentPatterns call completes quickly. This call only occurs once, and therefore optimising it will not significantly improve the speed of the code.

However, the outer for loop takes a great deal of time to complete (in the order of several hours). In particular, the patternOccurances call is particularly slow. However, I am sure other inefficiencies exist in the code.

Other programmers have completed the same task within 40s.

• How can my coding be improved?
• How can I optimise my algorithms performance?
• Do more efficient algorithms exist?

For those interested, the large dataset for which the above algorithms should efficiently work can be found here (4.2MB).

Answer 1

1. Introduction

This code has so many basic problems with documentation, layout, spelling, naming, and so on, that I didn't get around to looking at what it does. In a professional context we would say, "this code is not remotely ready for review" and send it back for you to fix!

It's worth reading the Python style guide (PEP8) and following the advice there. It's not compulsory, but it will make it easier for you to collaborate with other Python programmers and read each others' code. In particular, I'd use 4 spaces indentation, and lowercase_words_with_underscores for function and method names.

2. The function patternOccurances2

1. Your main code calls patternOccurances but when you define the function it's called patternOccurances2.

2. There's no documentation for patternOccurances2. What does it do and how am I supposed to call it?

3. This function takes default arguments text="GATATATGCATATACTT", pattern="ATAT". These do not seem like suitable default arguments. Default arguments are normally used for common values, to save the caller having to specify them. If you want to provide test data, then write a test case separately. See below for one way to do it.

4. This function takes its arguments in the order (text, pattern) but it would be better for it to take the arguments the other way round, for consistency with functions in the standard library like re.search.

5. This function starts by calling:

   position = config.genome.find(pattern, position+1)
   

   What is config.genome and how is it relevant here? It looks as if it might be a mistake for text.

6. It seems inelegant to append a -1 to the list of occurrences, and then remove it again. It would be better to write the loop like this:

   result = []
   position = -1
   while True:
       position = text.find(pattern, position + 1)
       if position == -1:
           return result
       result.append(position)
   

7. Instead of returning a list, it's nearly always better in Python to generate the results using the yield instruction, explained here. This means that you can process the results one by one without having to build up an intermediate list in memory. So I would write the function like this:

   def pattern_occurrences(pattern, text):
       """Generate the indices of the (possibly overlapping) occurrences of
       pattern in text. For example:
   
       >>> list(pattern_occurrences('ATAT', 'GATATATGCATATACTT'))
       [1, 3, 9]
   
       """
       position = -1
       while True:
           position = text.find(pattern, position + 1)
           if position == -1:
               return
           yield position
   

   Notice that the docstring includes an example test case. You can run this and check that the result is correct using the doctest module.

3. The function frequentPatterns

1. Something's gone wrong with the indentation in this function. A copy-paste error?

2. The documentation says Input ... number of occurances t. But the actual argument is called numberOfOccurances, not t.

3. The documentation for frequentPatterns says Output: All most frequent k-mers in text. But actually, only the k-mers that appear at least t times are output.

4. There's no explanation in the documentation for the behaviour if you don't specify a value for t.

5. The documentation for frequentPatterns is in a comment. But that means that I can't read it from the interactive interpreter:

   >>> help(frequentPatterns)
   Help on function frequentPatterns in module __main__:
   
   frequentPatterns(text='', k=2, numberOfOccurances=None)
   

   If you put this documentation into a docstring, like this:

   def frequent_patterns(text="", k=2, t=None):
       """Return a list of the most frequent k-mers in a string.
   
       Input: A string text, an integer k, and a number of occurrences t.
       Output: All k-mers that appear at least t times in text.
   
       If t is omitted, return a list containing the k-mer that appears
       the most often (or the k-mers, if there is a tie).
   
       """
       words = {}
       # ...
   

   then you'd be able to read it from the interactive interpreter:

   >>> help(frequent_patterns)
   Help on function frequent_patterns in module __main__:
   
   frequent_patterns(text='', k=2, t=None)
       Return a list of the most frequent k-mers in a string.
   
       Input: A string text, an integer k, and a number of occurrences t.
       Output: All k-mers that appear at least t times in text.
   
       If t is omitted, return a list containing the k-mer that appears
       the most often (or the k-mers, if there is a tie).
   

Answer 2

"Occurrances" should be spelled as "occurrences".

The core of your strategy is this loop from frequentPatterns():

for i in range(len(text)-k + 1): #iterate over the valid length of the string
  a = text[i:i+k]
  if a in words:
    words[a] = words[a] + 1
  else:
    words[a] = 1

I suggest using collections.Counter to keep track of the word count, as that is exactly what it is for.

words = collections.Counter([text[i:i+k] for i in range(len(text) - k + 1)])

Taking all those substrings might be a performance concern. Furthermore, you will have to find those substrings later in patternOccurances(). If k ≤ 16, you will likely improve performance by converting the words into numbers. Since there are four nucleobases, it should take only 2 bits to store each nucleobase. You should be able to store a 16-nucleobase string in a 32-bit number.

Here's my attempt at it:

from array import array
from collections import Counter

DNA_NUCLEOBASES = 'ACTG'

def numerify(s, k):
    '''
    Summarizes consecutive k-length substrings of s as as list of numbers.
    s is a string consisting of letters from DNA_NUCLEOBASES.
    k is an integer, 1 <= k <= 16.
    '''
    TRANS = dict(map(lambda c: (c, DNA_NUCLEOBASES.index(c)), DNA_NUCLEOBASES))
    mask = 4 ** k - 1
    # Depending on k, we might need a byte, short, or long.
    array_type = '!BBBBHHHHLLLLLLLL'[k]
    result = array(array_type, [0] * (len(s) - (k - 1)))
    v = 0
    for i in range(len(s)):
        v = (v << 2 & mask) | TRANS[s[i]]
        result[i - (k - 1)] = v
    return result.tolist()

def stringify(n, k):
    '''
    Converts a number from numerify() back into a k-length string of
    DNA_NUCLEOBASES.
    '''
    result = [DNA_NUCLEOBASES[(n >> 2 * i) & 3] for i in range(k)][::-1]
    return ''.join(result)

def find_clumps(s, k):
    '''
    Lists all k-length substrings of s in descending frequency.
    Returns a list of tuples of the substring and a list of positions at which
    they occur.
    '''
    clumps = numerify(s, k)
    result = []
    for clump, occurrences in Counter(clumps).most_common():
        if occurrences <= 1:
            break
        positions = []
        i = -1
        for _ in range(occurrences):
            i = clumps.index(clump, i + 1)
            positions.append(i)
        result.append((stringify(clump, k), positions))
    return result


s = 'CGGACTCGACAGATGTGAAGAACGACAATGTGAAGACTCGACACGACAGAGTGAAGAGAAGAGGAAACATTGTAA'
k = 5
for clump, positions in find_clumps(s, k):
    print("%s occurs at %s" % (clump, positions))

I'll leave it to you to implement the L-length windows.

Answer 3

Here is my version, which takes a little over 6 seconds to run in my tests.

The main goal of this version is to avoid iterating over the same piece of text repeatedly. In more naive implementations you might start with the first 'L' window of the genome, and check all of its contents before moving on to the next 'L' size window.

However, this window is only really changing by two characters - one is being removed from the front, and one added to the end. This means that we can keep the substring counts from the first window, decrement the counter for the first substring (as it's no longer there) and increment the counter for the new final string. This proves to be much faster that repeatedly iterating and counting.

Further improvements might be to discard the idea of storing the L size window, in case a future analysis involves a very large window where storing an extra copy of it would be memory intensive.

My initial version used a Counter, but this proved to be slower than a defaultdict, so I removed it. The list of relevant substrings is also more intuitively stored as a list, but the lookup of ... in list proved to scale poorly when the size of relevant matches grew large, so that is stored as a dictionary for better performance.

#!/usr/bin/env python

from collections import defaultdict

with open('E-coli.txt', 'r') as f:
    genome = f.read().rstrip()


kmer = 9  # Length of the k-mer
windowsize = 500  # genome substring length to register clumps in
min_clumpsize = 3  # minimum number of repetitions of the k-mer


# find substrings at least kmer long that occur at least min_clumpsize
# times in a window of windowsize in genome

def get_substrings(g, k):
    """
Take the input genome window 'g', and produce a list of unique 
substrings of length 'k' contained within it. 
    """
    substrings = list()

    # Start from first character, split into 'k' size chunks
    # Move along one character and repeat. No sense carrying on beyond
    # a starting point of 'k' since that will be the first iteration again.
    for i in range(k):
        line = g[i:]
        substrings += [line[i:i + k]
                       for i in range(0, len(line), k) if i + k <= len(line)]

    # Using collections.Counter increases the runtime by about 3 seconds,
    # during testing.
    results = defaultdict(int)
    for s in substrings:
        results[s] += 1
    return results


def find_clumps(genome, kmer, windowsize, clumpsize):
    """
In a given genome, examines each windowsize section for strings of length kmer
that occur at least clumpsize times. 

Input: 
genome: text string to search
kmer:  length of string to search for
windowsize: size of the genome section to consider for clumping
clumpsize: the kmer length strings must occur at least this many times

Returns: a list of the strings that clump
    """
    window = genome[0:windowsize]

    # Initialise our counter, because the main algorithm can't start from
    # scratch.
    patterns = get_substrings(window, kmer)

    # Using a dictionary not a list because the lookups are faster once the
    # size of the object becomes large
    relevant = {p: 1 for p in patterns if patterns[p] >= clumpsize}

    starting_string = genome[0:kmer]

    for i in range(windowsize, len(genome)):
        # Move the window along one character
        window = window[1:]
        window += genome[i]

        # This is the only string that can decrease if we've moved one
        # character
        patterns[starting_string] -= 1
        starting_string = window[0:kmer]

        # This is the only string that can increase if we've moved one
        # character
        ending_string = window[-kmer:]
        patterns[ending_string] += 1

        # if there are enough matches of the string at the end, add it to
        # matches.
        if patterns[ending_string] >= clumpsize and ending_string not in relevant:
            relevant[ending_string] = 1

    return list(relevant)


if __name__ == "__main__":
    clumps = find_clumps(genome, kmer, windowsize, min_clumpsize)
    print("Total: {}".format(len(clumps)))

Answer 4

Finding frequent occurring patterns in a string is a special case of the problem of finding all patterns that find more than a specified number of times in multiple strings. (Just set the specified number of times to your threshold for frequent and the number of strings to be searched to 1.)

Finding all patterns that find more than a specified number of times in multiple strings is addressed here. This link contains Python and C++ solutions of increasing complexity and speed that find all the repeated patterns that I need to find in the strings I work with. I think all these solutions can be easily adapted to your specific problem. Let me know if they can't

I have heard that bioinformaticians use compressed suffix trees to solve this sort of problem. You may find some faster solutions by following that line of investigation.