Finding Clumps in Python

Learning Objectives

In the previous code along, we used dictionaries to implement a frequency table in order to find the most frequent words of given length in a text. Now that we are comfortable with dictionaries, we will return to the problem of finding clumps in a genome, a problem that can help us scan the genome for areas that have surprisingly frequent k-mers, which may imply that these k-mers serve as “hidden messages” (i.e., that they are involved in protein-DNA binding).

Recall from the core text that we say that a k-mer pattern forms an (L, t)- clump inside a (longer) string genome if there is an interval of genome of length L in which this k-mer appears at least t times. We now restate the Clump Finding Problem.

We then introduced pseudocode for an algorithm called find_clumps(), which ranges over every substring window of length L of a given string text. For each such string window, it builds a frequency table for all substrings of length k. Any substrings in this table occurring t or more times are clumps, and we will add them to some list patterns if they do not already occur in patterns.

FindClumps(text, k, L, t)
    patterns ← an array of strings of length 0
    n ← length(text)
    for every integer i between 0 and n − L
        window ← text[i, i + L]
        freqMap ← FrequencyTable(window, k)
        for every key s in freqMap
            if freqMap[s] ≥ t and Contains(patterns, s) = false
                patterns ← append(patterns, s)
    return patterns

Contains(patterns, s)
    for every string pattern in patterns
        if s = pattern
            return true
    return false

In this lesson, we will implement FindClumps() along with the Contains() subroutine that identifies whether a given slice of strings patterns contains a string s. We will then apply FindClumps() to the E. coli genome to identify all short strings that occur frequently in short regions and may serve as hidden messages to the cell.

Code along summary

Setup

Create a folder called clumps in your python/src directory and create a text file called main.py in the python/src/clumps folder. We will edit main.py, which should have the following starter code.

def main():
    print("Clumps with Python!")

if __name__ == "__main__":
    main()

A first implementation of finding clumps

We begin with implementing find_clumps(). We start with some parameter checking. We will also change the window length L from the pseudocode into a more descriptive parameter window_length.

def find_clumps(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Output:
    - list: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

Because we don’t know the number of strings that will form clumps in advance, we first create a list of strings patterns of length zero, to which we will append any clumps that we find. We will eventually return this list patterns.

def find_clumps(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Output:
    - list: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # To fill in...

    return patterns

We next will fill in the meat of the find_clumps() function. To do so, we will need to range over every possible “window” substring of text having length window_length, forming the frequency table for each such window. We know from our work with ranging over all substrings of a text of given length that the starting positions of all these windows will range from 0 to n - window_length, inclusively, where n is the length of text. We can now fill in the beginning of our for loop.

def find_clumps(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Output:
    - list: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # range over all possible windows of text
    for i in range(0, n - window_length + 1):
        # set the current window
        window = text[i:i + window_length]

        # let's make the frequency table for this window
        freq_map = frequency_table(window, k)

    # To fill in...

    return patterns

We now will consider what to do once we have formed the frequency table freq_map of a given window. We want to range through this dictionary and grab all k-mer keys whose corresponding integer values are at least equal to t, the threshold value for considering a string to be a clump. If we find such a string, and it does not already occur in patterns, then we should append this string to patterns. We now have completed our find_clumps() function.

def find_clumps(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Output:
    - list: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # range over all possible windows of text
    for i in range(0, n - window_length + 1):
        # set the current window
        window = text[i:i + window_length]

        # let's make the frequency table for this window
        freq_map = frequency_table(window, k)

        # what occurs frequently (i.e., t or more times)?
        for s, val in freq_map.items():
            # append s to patterns if s occurs frequently and s doesn't already appear in patterns
            if val >= t and contains(patterns, s) == False:
                patterns.append(s)

    return patterns

As for contains(), this type of subroutine might start feeling comfortable now. To implement this function, we can use our in operator to find if a an element is contained in a list. If we ever find a value of the slice patterns that is equal to our query string s, then we return true, and if we range over all values of patterns without finding any value that is equal to s, then we return false.

However, we don’t need to write contains() at all, because in the previous code along, we introduced the in operator to check whether an element is a key within a dictionary. In this case, we can use the same operator to determine if an element is a value in a list (in this case, we are checking if it is not in the list). This allows us to update find_clumps() as follows, without using a contains() subroutine.

def find_clumps(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Output:
    - list: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # range over all possible windows of text
    for i in range(0, n - window_length + 1):
        # set the current window
        window = text[i:i + window_length]

        # let's make the frequency table for this window
        freq_map = frequency_table(window, k)

        # what occurs frequently (i.e., t or more times)?
        for s, val in freq_map.items():
            # append s to patterns if s occurs frequently and s doesn't already appear in patterns
            if val >= t and not (s in patterns):
                patterns.append(s)

    return patterns

Running the clump finding algorithm

We would like to apply find_clumps() to a real genome. First, let’s apply it to a smaller sample dataset just to make sure that everything appears to be in order. Using the following dataset, we should print only a single occurrence of "AA".

def main():
    print("Clumps with Python!")

    text = "AAAACGTCGAAAAA"
    k = 2
    window_length = 4
    t = 2

    print(find_clumps(text, k, window_length, t)) # Should print "AA".

Applying the clump finding algorithm to a bacterial genome, and installing modules

We are ready to scale up our work to a real bacterial genome to see how many strings appear surprisingly frequently in short regions of this genome. We will choose E. coli, the most commonly studied bacterium. You can view the E. coli genome as a .txt file at the Bioinformatics Algorithms website (click here to view). However, the genome has over 4.5 million nucleotides, and copying it into main() would pose a struggle. Instead, let’s write some code to read this genome from the URL over the internet.

First, however, we need to install the requests module, which will allow us to access a web page’s contents. To do so, we will use pip, the Python package installer, which is included in the installation of Python. To do so, open a command-line terminal and execute pip3 install requests (macOS/Linux) or pip install requests (Windows).

Then, include the following import at the top of main.py.

import requests

We next define the address that we wish to access from above as a string url, and we pass this string into a function from the requests module called get() that attempts to access the address. This function returns a “response” object (we will say more about objects about later in this course).

def main():
    # Code omitted for clarity

    url = "https://bioinformaticsalgorithms.com/data/realdatasets/Replication/E_coli.txt"

    response = requests.get(url)

We next ensure that there was no issue with opening the web page. To do so, we call a function raise_for_status() that will give an error if there was any problem. If not, we proceed to the next line, where we access response.text, which is a string representing the content of the web page; this string is exactly what we want, and so we declare a variable genome equal to this string.

def main():
    # Code omitted for clarity

    url = "https://bioinformaticsalgorithms.com/data/realdatasets/Replication/E_coli.txt"

    response = requests.get(url)

    response.raise_for_status() # give us an error if there was a problem 

    genome = response.text

We now will print the length of genome to ensure that it has been read in, which will show that it has over 4.6 million nucleotides. After we have done so, we establish some parameters, call find_clumps() on genome with these parameters, and print the number of resulting patterns that we found.

def main():
    # Code omitted for clarity

    url = "https://bioinformaticsalgorithms.com/data/realdatasets/Replication/E_coli.txt"

    response = requests.get(url)

    response.raise_for_status() # give us an error if there was a problem 

    genome = response.text
    
    print("The number of nucleotides in E. coli genome is", len(genome))
    
    k = 9
    window_length = 500
    t = 3

    clumps = find_clumps(genome, k, window_length, t)
    print("Found " + str(len(clumps)) + " total patterns as clumps.")

In our code, we found nearly 2,000 9-mers that are (500, 3)-clumps in the E. coli genome, meaning that a great deal of strings appear surprisingly frequently in short regions and may be involved as “hidden messages”. Yet before returning to the main text, we would like to observe that we can make find_clumps() faster.

Optimizing the clump finding algorithm by avoiding many subroutine calls

“Premature optimization is the root of all evil.”

Donald Knuth

In the light of Knuth’s famous quotation, we should always be careful pursuing optimizations before we have a working solution to our problems. Yet in this particular case, we already have an intuitive solution, and so our desire to optimize it is not premature.

The critical insight is that the frequency tables of consecutive substring “windows” in any string text are nearly identical. For example, say that we are working with the text "BANANASPLIT", the window length is equal to 6, and k is equal to 3. Then the frequency table of the first window ("BANANA") is shown in the table below.

Pattern	Count
"BAN"	1
"ANA"	2
"NAN"	1

The frequency table of the second window ("ANANAS") is shown in the table below. The only difference between this frequency table and the one above is that we lost one occurrence of "BAN", the first k-mer in "BANANA", and we gained one occurrence of "NAS", the final k-mer in "ANANAS".

Pattern	Count
"ANA"	2
"NAN"	1
"NAS"	1

In general, say that window_1 and window_2 are consecutive length-L substrings of a string text. The frequency table of window_1 will be the same as that of window_2, except that the former will include one additional occurrence of the first k-mer in window_1, and the latter will include one additional occurrence of the final k-mer in window_2.

To use this observation to optimize find_clumps(), we note that after generating the frequency table of the first window of text, we do not need to call frequency_table() to generate the frequency table of the second window of text. We need only to modify the frequency table of text[0:window_length] by decreasing the value associated with text[0:k] by 1, and increasing the value associated with text[window_length + 1 - k:window_length + 1] by 1.

More generally, after calling frequency_table() to form the frequency table of text[0:window_length], we will range an integer i between 1 and n - window_length + 1 (to access all remaining windows). For each such i, we update the current frequency table by decreasing the value in the frequency table associated with text[i - 1:i - 1 + k] and increasing the value associated with text[i + window_length - k:i + window_length].

This idea is implemented with the updated find_clumps_faster() function below. First, we perform some parameter checks and declare our list patterns (which we will eventually return) as well as our dictionary found_patterns.

def find_clumps_faster(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Returns:
    - list[str]: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # map to track whether I have added a string to patterns yet
    found_patterns = {}

    # to fill in

    return patterns

In the first window, every k-mer is being seen for the first time, so no check is needed as to whether we have encountered it; if a substring occurs t or more times within the window, then we add it to patterns.

def find_clumps_faster(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Returns:
    - list[str]: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # map to track whether I have added a string to patterns yet
    found_patterns = {}

    first_window = text[:window_length]
    freq_map = frequency_table(first_window, k)

    # append any patterns we find to patterns slice
    for s, freq in freq_map.items():
        if freq >= t:
            patterns.append(s)
            found_patterns[s] = True

    # range over all remaining possible windows of text
    # to fill in

    return patterns

We then range over all remaining windows. In each window, we first decrease by 1 the value associated with the first substring of length k in the preceding window, which is text[i - 1:i - 1 + k]. If the value associated with this string in freq_map becomes zero as a result, then we clean up the dictionary by using Python’s built-in del operator to delete this pattern from freq_map.

Then it is simply a matter of incrementing the frequency table associated with the substring at the end of the current window, which is text[i + window_length - k:i + window_length]. Once we have done so, we are ready to check any k-mers that occurred at least t times in the window, and add them to patterns if we have not yet seen them.

def find_clumps_faster(text: str, k: int, window_length: int, t: int) -> list[str]:
    """
    Finds a list of strings representing all k-mers that appear at least t times
    in a window of given length in the string.

    Parameters:
    - text (str): An input string.
    - k (int): k-mer's of size k.
    - window_length (int): the size of substrings of text in which we are looking for clumps
    - t (int): The k-mers must appear at least t amount of times.

    Returns:
    - list[str]: A list of k-mers that occur at least t times in a window of length window_length in text.

    """
    if len(text) == 0:
        raise ValueError("text must not be empty")
    if k <= 0 or window_length <= 0 or t <= 0:
        raise ValueError("k, window_length, and t must be positive")
    if k > len(text):
        return []
    if window_length < k:
        raise ValueError("window length must be at least equal to k")

    patterns: list[str] = []
    n = len(text)

    # map to track whether I have added a string to patterns yet
    found_patterns = {}

    first_window = text[:window_length]
    freq_map = frequency_table(first_window, k)

    # append any patterns we find to patterns slice
    for s, freq in freq_map.items():
        if freq >= t:
            patterns.append(s)
            found_patterns[s] = True

    # range over all remaining possible windows of text
    for i in range(1, n - window_length + 1):
        # decrease by 1 the value associated with the first substring of length k in the preceding window
        old_pattern = text[i - 1:i - 1 + k]
        freq_map[old_pattern] -= 1

        # clean up the map if the frequency of old_pattern was 1
        if freq_map[old_pattern] == 0:
            del freq_map[old_pattern]

        # add the pattern from the end of the current window
        new_pattern = text[i + window_length - k:i + window_length]
        freq_map[new_pattern] = freq_map.get(new_pattern, 0) + 1

        # I have updated the frequency map :)
        for s, freq in freq_map.items():
            is_key = found_patterns.get(s, False)
            if freq >= t and is_key == False:
                patterns.append(s)
                found_patterns[s] = True

    return patterns

Timing the optimized approach

Now that we have optimized our clump finding algorithm, let’s see if it truly is any faster by timing both find_clumps() and find_clumps_faster(). On our machine, the former runs in about 5 minutes, whereas the latter completes in under 50 seconds. Once again, we encounter the paradigm that we introduced in Chapter 0 that efficient algorithms are particularly important when applying them to large datasets.

def main():

    # code for reading in genome from file ...
    # (assume genome, k, window_length, t are defined above)

    start = time.time()
    _ = find_clumps(genome, k, window_length, t)
    elapsed1 = time.time() - start
    print(f"find_clumps {elapsed1:.6f} seconds")

    start = time.time()
    clumps = find_clumps_faster(genome, k, window_length, t)
    elapsed2 = time.time() - start
    print(f"find_clumps_faster {elapsed2:.6f} seconds")

    print(f"Speedup: {elapsed1/elapsed2:.2f}x faster")

Looking ahead

Unfortunately, although our clump finding algorithm finds candidate frequent words across the entire genome, it finds too many clumps, which means that if we are looking for a specific region like the origin of replication, we will need a targeted approach. We will introduce just such a method in the next code along.

Check your work from the code along

We provide autograders in the window below (or via a direct link) allowing you to check your work for the following functions:

• find_clumps()

powered by Advanced iFrame

---