Sequence Analysis with Python¶

The following assignments introduce applications of hashing with dict() primitive of Python. While doing so, a rudimentary introduction to biological sequences is given. This framework is then enhanced with probabilities, leading to routines to generate random sequences under some constraints, including a general concept of Markov-chains. All these components illustrate the usage of dict(), but at the same time introduce some other computational routines to efficiently deal with probabilities.
The function collections.defaultdict can be useful.

Below are some "suggested" imports. Feel free to use and modify these, or not. Generally it's good practice to keep most or all imports in one place. Typically very close to the start of notebooks.

The automated TMC tests do not test cell outputs. These are intended to be evaluated in the peer reviews. So it is still be a good idea to make the outputs as clear and informative as possible.

To keep TMC tests running as well as possible it is recommended to keep global variable assignments in the notebook to a minimum to avoid potential name clashes and confusion. Additionally you should keep all actual code exection in main guards to keep the test running smoothly. If you run check_sequence.py in the part07-e01_sequence_analysis folder, the script should finish very quickly and optimally produce no output.

If you download data from the internet during execution (codon usage table), the parts where downloading is done should not work if you decide to submit to the tmc server. Local tests should work fine.

DNA and RNA¶

A DNA molecule consist, in principle, of a chain of smaller molecules. These smaller molecules have some common basic components (bases) that repeat. For our purposes it is sufficient to know that these bases are nucleotides adenine, cytosine, guanine, and thymine with abbreviations A, C, G, and T. Given a DNA sequence e.g. ACGATGAGGCTCAT, one can reverse engineer (with negligible loss of information) the corresponding DNA molecule.

Parts of a DNA molecule can transcribe into an RNA molecule. In this process, thymine gets replaced by uracil (U).

1. Write a function dna_to_rna to convert a given DNA sequence $s$ into an RNA sequence. For the sake of exercise, use dict() to store the symbol to symbol encoding rules. Create a program to test your function.

Idea of solution¶

The solution involves converting a DNA sequence into an RNA sequence by replacing each occurrence of the nucleotide thymine (T) with uracil (U) and leaving all other nucleotides unchanged. This transformation is achieved using a dictionary that maps each nucleotide to its corresponding RNA counterpart. The function iterates over the input string, performs the replacement using the dictionary, and constructs the RNA sequence using Python's join method.

Discussion¶

The problem of converting a DNA sequence to an RNA sequence is relatively straightforward and can be efficiently handled using a dictionary for nucleotide mapping. Further improvements might include f.ex. input validation, but that doesn't appear necessary in this context.

Proteins¶

Like DNA and RNA, protein molecule can be interpreted as a chain of smaller molecules, where the bases are now amino acids. RNA molecule may translate into a protein molecule, but instead of base by base, three bases of RNA correspond to one base of protein. That is, RNA sequence is read triplet (called codon) at a time.

2. Consider the codon to amino acid conversion table in http://htmlpreview.github.io/?https://github.com/csmastersUH/data_analysis_with_python_2020/blob/master/Codon%20usage%20table.html. Write a function get_dict to read the table into a dict(), such that for each RNA sequence of length 3, say $\texttt{AGU}$, the hash table stores the conversion rule to the corresponding amino acid. You may store the html page to your local src directory, and parse that file.

Idea of solution¶

The get_dict() function parses a file containing codon-to-amino acid mappings and creates a dictionary from this data. The function uses a regular expression to find and extract codon-amino acid pairs from each line of the file. The regex pattern \b([AUGC]{3})\s+([A-Z*])\b captures RNA triplets (codons) and their corresponding amino acids. The function then constructs a dictionary with codons as keys and amino acids as values, and returns this dictionary.

Discussion¶

The initial approach of using simple string operations to extract codon-amino acid pairs was prone to errors due to inconsistent spacing and formatting in the data. By employing regular expressions, the revised solution provides a more robust and precise method for parsing the file. Debugging output helps verify that the dictionary is correctly populated. This method appears both efficient and adaptable, making it a reliable solution for the task.

3. Use the same conversion table as above, but now write function get_dict_list to read the table into a dict(), such that for each amino acid the hash table stores the list of codons encoding it.

Idea of solution¶

The function get_dict_list() aims to create a dictionary, mapping each amino to a list of codons that encode it.

• The process begins by reading a file containing codon-to-amino mappings.
• These mappings are extracted using the same regular expression pattern used in the last task.
• Each codon-amino pair is stored in the codon_to_amino_dict dictionary.
• Subsequently, the function converts this dictionary into the desired format, amino_to_codon_dict_list, which maps each amino to a list of codons.
• The conversion is done by iterating through the codon_to_amino_dict and appending codons to the appropriate list in amino_to_codon_dict_list.

The debug output, if enabled, provides intermediate verification of these mappings.

Discussion¶

The solution achieves the task by first populating a codon_to_amino_dict with codon-amino acid pairs and then transforming this dictionary into the final format where amino acids are keys and lists of codons are values. This two-step approach ensures clarity in the process but may involve some redundant steps, combining these steps into a single pass without intermediate storage could improve efficiency.

With the conversion tables at hand, the following should be trivial to solve.

4. Fill in function rna_to_prot in the stub solution to convert a given DNA sequence $s$ into a protein sequence. You may use the dictionaries from exercises 2 and 3. You can test your program with ATGATATCATCGACGATGTAG.

Idea of solution¶

The function rna_to_prot converts an RNA sequence into a protein sequence leveraging the codon-to-amino dictionary.

• First the DNA string is converted to RNA using earlier dna_to_rna function, and rna_to_prot function called with the result.
• In this function the codon-to-amino dictionary is fetched.
• Then the RNA sequence is split into codons (3-letter segments).
• Each codon is translated into its corresponding amino acid using the fetched dictionary.
  • This is achieved through list comprehension, which maps each codon to its corresponding amino and collects the aminos into a list.
• Finally, the list of amino acids is concatenated into a single string to form the protein sequence.

This approach ensures that each codon is accurately translated and the resulting protein sequence is produced.

Discussion¶

The rna_to_prot function translates an RNA sequence into a protein sequence in a straightforward manner by using the previously created functions. The use of list comprehension for splitting the RNA sequence into codons and then translating each codon into an amino acid provides a compact and readable solution. Debugging output helps verify the correctness of intermediate steps, by showing the list of codons and the resulting amino acids.

You may notice that there are $4^3=64$ different codons, but only 20 amino acids. That is, some triplets encode the same amino acid.

Reverse translation¶

It has been observed that among the codons coding the same amino acid, some are more frequent than others. These frequencies can be converted to probabilities. E.g. consider codons AUU, AUC, and AUA that code for amino acid isoleucine. If they are observed, say, 36, 47, 17 times, respectively, to code isoleucine in a dataset, the probability that a random such event is AUU $\to$ isoleucine is 36/100.

This phenomenon is called codon adaptation, and for our purposes it works as a good introduction to generation of random sequences under constraints.

5. Consider the codon adaptation frequencies in http://htmlpreview.github.io/?https://github.com/csmastersUH/data_analysis_with_python_2020/blob/master/Codon%20usage%20table.html and read them into a dict(), such that for each RNA sequence of length 3, say AGU, the hash table stores the probability of that codon among codons encoding the same amino acid. Put your solution in the get_probabability_dict function. Use the column "([number])" to estimate the probabilities, as the two preceding columns contain truncated values.

Idea of solution¶

The goal of this solution is to calculate the probability distribution of codons encoding each amino based on their frequency in a provided dataset. The program reads codon frequencies from a file and processes these frequencies to produce a normalized probability distribution for each codon-amino pairing.

The process consists of several steps:

• Read Codon Frequency Data: The program reads the codon frequency data from a file. This file contains lines with codons, amino acids, and their respective frequencies.

• Extract Frequencies Using Regular Expressions: A regular expression is used to parse each line of the file and extract codons along with their frequencies. The regular expression pattern matches a codon and its frequency value.

• Map Codons to Amino Acids: The program maps each codon to its corresponding amino, fetching the amino (using the previously defined dictionary, amino_to_codon_dict_list).

• Calculate Frequency Dictionary: For each codon, the program calculates the frequency for each aminos codon. These frequencies are stored in a nested dictionary structure (codon_adaptation_freq_dict), where the outer keys are aminos, and the inner dictionaries map codons to their frequencies.

• Compute Probabilities: Once all codon frequencies are collected, the program calculates the probability of each codon for each amino acid. This is done by dividing the frequency of each codon by the total frequency of all codons encoding the same amino. The result is a probability distribution stored in a new dictionary (codon_encoding_amino_probability_dict).

• Flatten the Probability Dictionary: The nested probability dictionary is then flattened for use in the main block of the program, resulting in a single dictionary that maps each codon directly to its probability.

• Output the Results: Finally, the program sorts the codon-probability pairs alphabetically by codon and prints them in a formatted table.

Discussion¶

The solution leverages a structured approach to convert raw frequency data into a usable probability distribution format. The use of regular expressions allows for flexible and robust parsing of the input data, ensuring that only correctly formatted lines are processed. The mapping from codons to aminos is managed using dictionaries, allowing quick lookups and updates.

A significant challenge encountered during development was managing inconsistent spacing and avoiding unnecessary data capture in the regex pattern. Initially, the regular expression was too broad, capturing more information than needed and struggling to handle varying spaces between data elements. To resolve this, the regex pattern was refined to focus solely on essential components: the three-letter codon sequence and the frequency count inside parentheses.

The program includes optional debug output. When enabled, the program prints detailed information about the internal state at various steps.

Now you should have everything in place to easily solve the following.

6. Write a class ProteinToMaxRNA with a convert method which converts a protein sequence into the most likely RNA sequence to be the source of this protein. Run your program with LTPIQNRA.

Idea of solution¶

The solution involves creating a class, ProteinToMaxRNA, that translates a protein sequence into its most likely corresponding RNA sequence based on codon usage probabilities. This is done using two key dictionaries: one that maps each amino to its possible codons (amino_to_codon), and another that provides the probability of each codon being used (codon_probabilities).

The convert method takes a protein sequence as input and processes each amino one by one. For each amino, the method retrieves its associated codons and their probabilities from the respective dictionaries. It then determines which codon has the highest probability of encoding the amino and appends this codon to a list. After processing all aminos in the input sequence, the method returns a concatenated string of the selected codons, forming the RNA sequence that is most likely.

Discussion¶

While the current implementation maintains a clear structure there may be some redundancy within the loop iteration. With some refinements, it could possibly be made even more efficient and concise. However, the code correctly implements the required functionality and provides the correct output for the example input and tests.

Now we are almost ready to produce random RNA sequences that code a given protein sequence. For this, we need a subroutine to sample from a probability distribution. Consider our earlier example of probabilities 36/100, 47/100, and 17/100 for AUU, AUC, and AUA, respectively. Let us assume we have a random number generator random() that returns a random number from interval $[0,1)$. We may then partition the unit interval according to cumulative probabilities to $[0,36/100), [36/100,83/100), [83/100,1)$, respectively. Depending which interval the number random() hits, we select the codon accordingly.

7. Write a function random_event that chooses a random event, given a probability distribution (set of events whose probabilities sum to 1). You can use function random.uniform to produce values uniformly at random from the range $[0,1)$. The distribution should be given to your function as a dictionary from events to their probabilities.

Idea of solution¶

The solution aims to randomly select an event from a given distribution where each event has a specified probability.

• First a random fraction between 0 and 1 is generated. This fraction will help determine which event to select based on its position relative to cumulative probabilities.

• The cumulative probabilities of the events are calculated from the given distribution. For this purpose, the list of individual event probabilities is transformed into a list where each element represents the sum of probabilities up to and including that event.

• Events and their corresponding cumulative values are paired together. This allows us to easily determine which event corresponds to the range that the random fraction falls within.

• By comparing the random fraction with the cumulative values, the correct event is selected. The first cumulative value that is greater than the random fraction indicates the event that should be returned.

If debugging is enabled, additional information is printed to help verify the correctness of the sampling process.

Discussion¶

The approach uses cumulative probabilities to map a random fraction to a specific event. This ensures that each event is selected in proportion to its specified probability. Debug output is integrated to aid in the verification of the sampling process. It provides visibility into key variables such as the distribution, random fraction, cumulative values, and the selected event.

The code works well with the given example input and the tests. For more genral use cases, there could be issues because the code assumes that the input probabilities are perfectly normalized to sum to exactly 1.0. One issue that came up while writing the code was list comprehension resulting in floating-point precision error, where cumulative values would come up such as 0.4499999999999999, instead of 0.45.

With this general routine, the following should be easy to solve.

8. Write a class ProteinToRandomRNA to produce a random RNA sequence encoding the input protein sequence according to the input codon adaptation probabilities. The actual conversion is done through the convert method. Run your program with LTPIQNRA.

Idea of solution¶

The solution involves mapping each amino in the input protein sequence to its possible codons using a dictionary (amino_to_codon). Another dictionary (codon_probabilities) provides the probabilities of each codon being used. The program iterates through each amino in the protein sequence, retrieves the possible codons and their associated probabilities, and randomly selects a codon using the random_event function. The selected codons are concatenated to form the final RNA sequence.

Debugging statements are included to trace the steps and verify correct function at each stage.

Discussion¶

This task was a relatively straightforward implementation of the previously created functions and the class in the last exercise. The print statements from the last exercise attempt to further visualize that when calling the random event class the right codon is returned, based on which fraction variable the random event class produced.

The task only required a straightforward implementation, using pre-existing dictionaries for amino mapping and codon probability data. The main challenge was ensuring that each amino acid was correctly converted into a corresponding RNA codon based on its probability distribution. The debug output effectively shows the mapping and selection process, confirming that the random_event function correctly chooses codons according to the provided probabilities.

Future improvements could include optimizing the handling of input data and refining the random selection process to ensure performance across larger datasets.

Generating DNA sequences with higher-order Markov chains¶

We will now reuse the machinery derived above in a related context. We go back to DNA sequences, and consider some easy statistics that can be used to characterize the sequences. First, just the frequencies of bases $\texttt{A}$, $\texttt{C}$, $\texttt{G}$, $\texttt{T}$ may reveal the species from which the input DNA originates; each species has a different base composition that has been formed during evolution. More interestingly, the areas where DNA to RNA transcription takes place (coding region) have an excess of $\texttt{C}$ and $\texttt{G}$ over $\texttt{A}$ and $\texttt{T}$. To detect such areas a common routine is to just use a sliding window of fixed size, say $k$, and compute for each window position $T[i..i+k-1]$ the base frequencies, where $T[1..n]$ is the input DNA sequence. When sliding the window from $T[i..i+k-1]$ to $T[i+1..i+k]$ frequency $f(T[i])$ gets decreases by one and $f(T[i+k])$ gets increased by one.

9. Write a generator sliding_window to compute sliding window base frequencies so that each moving of the window takes constant time. We saw in the beginning of the course one way how to create generators using generator expression. Here we use a different way. For the function sliding_window to be a generator, it must have at least one yield expression, see https://docs.python.org/3/reference/expressions.html#yieldexpr.

Here is an example of a generator expression that works similarily to the built in range generator:

def range(a, b=None, c=1):
    current = 0 if b == None else a
    end = a if b == None else b
    while current < end:
        yield current
        current += c

A yield expression can be used to return a value and temporarily return from the function.

Idea of solution¶

The sliding_window generator function calculates nucleotide frequencies within a sliding window of a specified size over a DNA sequence. The function first initializes a frequency dictionary for the initial window, ensuring all possible nucleotides ('A', 'C', 'G', 'T') are accounted for, even those with zero counts. The initial frequency dictionary is then yielded. As the window slides one position at a time through the sequence, the function updates the counts by decrementing the frequency of the nucleotide that exits the window and incrementing the frequency for the nucleotide that enters the window. This update process is efficient, operating in constant time since only two nucleotide frequencies change per step.

Discussion¶

Initially, the tests failed because the frequency dictionaries did not account for all nucleotide bases ('A', 'C', 'G', 'T'), specifically those that had zero frequencies. This was corrected by initializing the frequency dictionary with all nucleotides. The concept of sliding a window and adjusting counts in constant time was somewhat challenging to grasp at first. I initially included loops that would have caused the algorithm to operate in inconsistent time, which would have been inefficient. After understanding how to adjust only the counts of the nucleotides entering and exiting the window, the function now correctly maintains constant time for each window move.

The debug print statements are included to visually demonstrate the changes in nucleotide frequencies as the window slides across the DNA sequence. This helped confirm the correctness of the solution and provided insight into the algorithm's step-by-step operations.

Our models so far have been so-called zero-order models, as each event has been independent of other events. With sequences, the dependencies of events are naturally encoded by their contexts. Considering that a sequence is produced from left-to-right, a first-order context for $T[i]$ is $T[i-1]$, that is, the immediately preceding symbol. First-order Markov chain is a sequence produced by generating $c=T[i]$ with the probability of event of seeing symbol $c$ after previously generated symbol $a=T[i-1]$. The first symbol of the chain is sampled according to the zero-order model.
The first-order model can naturally be extended to contexts of length $k$, with $T[i]$ depending on $T[i-k..i-1]$. Then the first $k$ symbols of the chain are sampled according to the zero-order model. The following assignments develop the routines to work with the higher-order Markov chains. In what follows, a $k$-mer is a substring $T[i..i+k-1]$ of the sequence at an arbitrary position.

10. Write function context_list that given an input DNA sequence $T$ associates to each $k$-mer $W$ the concatenation of all symbols $c$ that appear after context $W$ in $T$, that is, $T[i..i+k]=Wc$. For example, GA is associated to TCT in $T$=ATGATATCATCGACGATGTAG, when $k=2$.

Idea of solution¶

The solution uses the sliding window approach to examine each k-mer (substring of length k) within the given sequence s.

• For each k-mer, it captures the k-mer and the character that immediately follows it in the sequence.
• These k-mers serve as keys in a dictionary (context_dict), where each key's value is a concatenated string of all characters that appeared after this k-mer throughout the sequence.
• If a k-mer already exists in the dictionary, the following character is concatenated to its existing values.
• Otherwise, a new entry is created with the k-mer as the key and the subsequent character as its value.

This approach efficiently builds a dictionary of k-mer contexts using simple string slicing and dictionary operations.

Discussion¶

The implementation of the context_list function is straightforward and effectively solves the problem by utilizing a simple loop to manage the sliding window across the string s. By directly manipulating the dictionary to store or update the context of each k-mer, the code avoids unnecessary complexity and remains readable.

Debug output allows for step-by-step tracing of the program’s operation.

The solution is again quite barebones, in that it doesn't check extensively for edge-cases or errors. For example, the function does not explicitly handle cases where the length of the string s is less than k. In such cases, the function would return an empty dictionary since there are no valid k-mers to process. Although it is the correct behaviour in this instance, it is an example of how edge cases could affect the programs process.

11. With the above solution, write function context_probabilities to count the frequencies of symbols in each context and convert these frequencies into probabilities. Run context_probabilities with $T=$ ATGATATCATCGACGATGTAG and $k$ values 0 and 2.

Idea of solution¶

The solution creates a dictionary of k-mers (substrings of length k) and the probability of subsequent characters based on a given sequence. The function context_probabilities:

• Uses function context_list to generate a dict that maps each k-mer in the sequence s to the characters that immediately follow it in the sequence
• For each k-mer, it counts the frequency of each character that follows it, storing these counts in a separate frequency dictionary (freq_dict).
• These frequencies are then converted into probabilities by dividing each frequency by the total number of occurrences for that k-mer.
• The result is stored in context_probability_dict, where each k-mer is a key, and its value is another dictionary containing characters and their corresponding probabilities.

The function returns a dictionary which provides a probabilistic model of the sequence's structure based on k-mer contexts.

Discussion¶

This exercise extends the logic of the previous exercise by calculating the probabilities of each character following specific k-mers. The implementation is straightforward and reuses previous logic with slight modifications. The debug output matches the expected results, demonstrating that the function correctly calculates probabilities. Further enhancements could include handling more edge cases, such as when k is larger than the length of s.

12. With the above solution and the function random_event from the earlier exercise, write class MarkovChain. Its generate method should generate a random DNA sequence following the original $k$-th order Markov chain probabilities.

Idea of solution¶

The solution has a class MarkovChain that generates a random DNA sequence based on the given Markov chain probabilities.

• The class is initialized with a zeroth probability dictionary (zeroth) and a k-th order transition probability dictionary (kth).

• The generate method is the core of the programme. It starts by optionally setting a random seed for reproducibility, ensuring that the sequence generation can be replicated.

• The first k characters of the sequence are generated using the zeroth-order probabilities, as the Markov chain cannot determine transitions without an initial state.

• After initializing the first k characters, the method generates the remaining characters using k-th order probabilities.

• It does this by looking at the last k characters (the current k-mer) and choosing the next character based on the transition probabilities defined in kth.

• If a k-mer does not exist in kth, the algorithm falls back to the zeroth-order probabilities, ensuring that the sequence generation continues.

Discussion¶

Initially, the problem seemed straightforward, but the tests made it more challenging. The requirement to handle missing k-mers in the kth dictionary was not immediately obvious from the task description. The approach taken uses a while loop to manage sequence generation which ensures flexibility. The use of a fallback mechanism for missing k-mers makes the solution robust against incomplete transition matrices.

Although not strictly related to the functioning of the solution, a potential improvement could be further conciseness of debug statements. These make the code more verbose, but i find the output helpful to visualize the sequence of events.

If you have survived so far without problems, please run your program a few more times with different inputs. At some point you should get a lookup error in your hash-table! The reason for this is not your code, but the way we defined the model: Some $k$-mers may not be among the training data (input sequence $T$), but such can be generated as the first $k$-mer that is generated using the zero-order model.

A general approach to fixing such issues with incomplete training data is to use pseudo counts. That is, all imaginable events are initialized to frequency count 1.

13. Write a new solution context_pseudo_probabilities based on the solution to problem 11. But this time use pseudo counts in order to obtain a $k$-th order Markov chain that can assign a probability for any DNA sequence. You may use the standard library function itertools.product to iterate over all $k$-mer of given length (product("ACGT", repeat=k)).

Idea of solution¶

The task requires implementing a pseudo-count based, context probability model, to handle unseen k-mers in a DNA sequence. The code has different approches based on the value of k:

• k = 0:

  • Create a dictionary for nucleotide frequencies, initialized with a pseudo-count.
  • Increment the counts based on the given sequence.
  • Compute probabilities by normalizing counts.

• k > 0:

  • Generate all possible k-mers using itertools.product.
  • Initialize k-mers with pseudo-counts.
  • Update k-mer counts based on the given sequence by sliding through the string.
  • Compute probabilities for each k-mer by normalizing counts.

This approach ensures that all possible k-mers have a non-zero probability, even if they were not observed in the training data.

Discussion¶

Initially, the task was difficult for me to complete, as I wasn't sure how to proceed or to satisfy all the tests. I tried using earlier functions to provide a base for completing this task, but all such implementations quickly became overly complex.

The solution addresses the issue of unseen k-mers by initializing all possible k-mers with a pseudo-count of 1. This approach prevents zero probabilities for k-mers not present in the sequence. The results from the provided example demonstrate that the function correctly computes both zeroth and kth order probabilities. Further testing and handling for edge cases could probably still improve the code.

14. Write class MarkovProb that given the $k$-th order Markov chain developed above to the constructor, its method probability computes the probability of a given input DNA sequence.

Idea of solution¶

The solution uses class MarkovProb to calculates the probability of a DNA sequence based on a given Markov model. Again, The code has different approches based on the value of k:

• For k = 0: The method iterates through each nucleotide in the sequence, fetches its probability from the zeroth-order probabilities andding it to a list, and computes the product of these probabilities.

• For k > 0: The method starts by calculating the probability of the initial k-mer based on individual nucleotide probabilities. It multiplies this with the overall probabiity value, then it iteratively calculates the transition probabilities for each subsequent nucleotide following the k-mers. The total probability is the product of the initial k-mer probability and all transition probabilities.

Discussion¶

During implementation, an issue was encountered for cases where k > 0.

Initially, in cases where k > 0, the computation didn't correctly account for the initial k-mer's probability, leading to incorrect results. This was corrected by separately calculating the probability for the initial k-mer and then using it as a base to add with subsequent transitions.

The use of functional constructs helped make the code more concise and readable. However, this also introduces a potential risk for debugging since these constructs can obscure the detailed flow of calculations. I tried to find a balance between these in the code.

Overall, while the code works correctly for the given inputs and tests while also being relatively optimized. More comprehensive testing across various sequences, k-values and probabilities could probably increase robustness.

With the last assignment you might end up in trouble with precision, as multiplying many small probabilities gives a really small number in the end. There is an easy fix by using so-called log-transform. Consider computation of $P=s_1 s_2 \cdots s_n$, where $0\leq s_i\leq 1$ for each $i$. Taking logarithm in base 2 from both sides gives $\log _2 P= \log _2 (s_1 s_2 \cdots s_n)=\log_2 s_1 + \log_2 s_2 + \cdots \log s_n= \sum_{i=1}^n \log s_i$, with repeated application of the property that the logarithm of a multiplication of two numbers is the sum of logarithms of the two numbers taken separately. The results is abbreviated as log-probability.

15. Write class MarkovLog that given the $k$-th order Markov chain developed above to the constructor, its method log_probability computes the log-probability of a given input DNA sequence. Run your program with $T=$ ATGATATCATCGACGATGTAG and $k=2$.

Idea of solution¶

The goal of this solution is to calculate the log-probability of a given DNA sequence based on the Markov model. Here's how the solution works:

• Initialize the Log Probability: Start with a log-probability of 0.0, which is equivalent to a probability of 1 in log space.

• Compute Initial k-mer Log-Probability: The probability of observing the initial k-mer (the first k characters of the sequence) is computed by summing the log2 probabilities of each character's occurrence, based on the zeroth-order probabilities.

• Iteratively Calculate Transition Log-Probabilities: For each subsequent nucleotide in the sequence, compute the conditional log-probability given the preceding k-mer using the k-th order probabilities. Sum these log-probabilities to the total log-probability.

• Combine Results: The total log-probability is obtained by summing the initial k-mer log-probability with the sum of transition log-probabilities. The final log-probability is rounded and returned as the result.

Discussion¶

The program follows a similar pattern as the last task, now calculating the log probability. It performs well with the given example sequence ATGATATCATCGACGATGTAG and k=2. The calculated log-probability matches expected values based on the provided zeroth and k-th order probabilities. It also passes all the tests.

The main problem I had with this task was understanding that I should have used log2 to get the correct results, which led quite alot of debugging, although it was clear from reading the task again. A difference compared to the previous task is that there isn't any code sperately handling input where k is 0. It didn't appear to be a requirement from the task description and the tests passed without it, so to keep the code more concise I left that part out.

Finally, if you try to use the code so far for very large inputs, you might observe that the concatenation of symbols following a context occupy considerable amount of space. This is unnecessary, as we only need the frequencies.

16. Optimize the space requirement of your code from exercise 13 for the $k$-th order Markov chain by replacing the concatenations by direct computations of the frequencies. Implement this as the better_context_probabilities function.

Idea of solution¶

The test already passed with my original code from task 13. In that implementation - context_dict = context_list(s, k) - concatenating the symbols to string denoting the full context is not called (as it is in exercise 11). In light of this, the solution contains some minor improvements, marked by the comment space improvement above, and detailed further below.

• ORIGINAL:

  • counts_dict = {nucleotide: pseudo_count for nucleotide in nucleotides}

• OPTIMIZED:

  • counts_dict = defaultdict(lambda: pseudo_count)

• BENEFIT:

  • By using a defaultdict with a default value (pseudo_count), we avoid pre-initializing each nucleotide key. This approach automatically initializes any accessed key with pseudo_count, thus saving memory that would be otherwise used for storing all potential keys upfronT.

• ORIGINAL:

  • kmers = [''.join(kmer) for kmer in product(nucleotides, repeat=k)]
  • counts_dict = {kmer: {nucleotide: pseudo_count for nucleotide in nucleotides} for kmer in kmers}

• OPTIMIZED:

  • counts = {''.join(kmer): {n: pseudo_count for n in nucleotides} for kmer in product(nucleotides, repeat=k)}

• BENEFIT:

  • The optimized approach directly initializes the count dictionary for all possible $k$-mers. By using this method, we avoid creating an excessive number of intermediate objects and maintain a leaner data structure that only contains the essential keys and values.

Discussion¶

As mentioned in the description of the solution I wasn't quite sure how to proceed with this task, as the original implementation already passed the tests. I think there would have been a wider range of different solutions I might have considered if it wasn't necessary to initialize pseudo-counts as well. As the task and the tests were set out however I felt a bit restricted in my ability to make significant changes, and settled for making a few minor improvements. The changes ensure that the solution effectively optimizes memory usage for both k == 0 and k > 0 cases by directly computing frequencies instead of storing extensive intermediate data.

While the earlier approach of explicit concatenation of symbols following a context suffered from inefficient use of space, it does have a benefit of giving another much simpler strategy to sample from the distribution: observe that an element of the concatenation taken uniformly randomly is sampled exactly with the correct probability.

17. Revisit the solution 12 and modify it to directly sample from the concatenation of symbols following a context. The function np.random.choice may be convenient here. Implement the modified version as the new SimpleMarkovChain class.

Idea of solution¶

The solution to this task is based on directly using the concatenation of symbols following a given context to sample from the distribution efficiently.

• It begins by initializing the random seed using np.random.seed(seed) to ensure reproducibility.

• For the special case when k = 0, the context length is zero, which means there are no preceding characters to consider. In this case, the solution generates n random characters directly from the input string s. This is achieved by uniformly sampling characters using np.random.choice.

• When k > 0, the solution starts by generating an initial string of k random characters from the input.

• It then constructs a list of all possible k+1-length substrings (transition_list) from the input string s. These substrings represent all possible transitions, where each transition is a k-mer followed by a single character that can follow that k-mer.

• The core of the solution is a loop that extends the generated sequence character by character.

  • At each step, the last k characters of the generated string are used to form the current k-mer.
  • The solution then finds all possible next characters from the transition_list that could follow this k-mer.
  • If possible transitions exist, the next character is sampled randomly from this list.
  • If no transitions are found, a fallback character is randomly chosen from the initial list of all characters in s.

By directly sampling from the list of possible transitions rather than maintaining a complex probability distribution or dictionary, the solution efficiently generates sequences according to the specified Markov chain order.

Discussion¶

The solution correctly handles both the case where k = 0 and cases where k > 0, generating sequences based on the input string's statistical properties. For k = 0, the output is a random sequence of characters from s, which matches expectations since there is no dependency on previous characters.

For k > 0, the generated sequence reflects the input string's character transitions for the specified k-mer length. The debug output shows the step-by-step generation process, including the construction of the transition list and the selection of each character. This helps ensure the logic aligns with the Markov chain principles.

During development, one of the challenges was handling cases where no valid transitions were found for a given k-mer, which could lead to an empty transition list and errors. The current solution mitigates this by including a fallback mechanism that samples from the original input string, thus ensuring that the generation process continues regardless of the transition list content.

$k$-mer index¶

Our $k$-th order Markov chain can now be modified to a handy index structure called $k$-mer index. This index structure associates to each $k$-mer its list of occurrence positions in DNA sequence $T$. Given a query $k$-mer $W$, one can thus easily list all positions $i$ with $T[i..k-1]=W$.

18. Implement function kmer_index inspired by your earlier code for the $k$-th order Markov chain. Test your program with ATGATATCATCGACGATGTAG and $k=2$.

Idea of solution¶

The objective of this task is to construct a k-mer index for a given DNA sequence s. The function kmer_index is implemented to achieve this:

• Handle Special Case (k = 0): The function returns a dictionary with an empty string as the key and a list containing all possible positions in the sequence. This is because a 0-mer effectively means "no context" so every position is a valid occurrence.

• General Case (k > 0): For values of k greater than 0, the function iterates over the sequence s to extract all possible k-mers. This is done by sliding a window of length k across the string:

  • For each position i from 0 to len(s) - k + 1, a k-mer starting at that position is extracted.
  • This k-mer is then added to a dictionary, kmer_indices_dict, with its starting index i appended to the list of indices for that k-mer.

Debugging output is provided to print the final k-mer index dictionary and its contents, which helps in verifying correctness and understanding the distribution of k-mers in the sequence.

Discussion¶

The solution creates a k-mer index by directly mapping each k-mer to its list of starting positions in the input string s.

Initially, the solution didn't account for k = 0, resulting in test failures and incorrect outputs. By explicitly handling this case, where every position in the sequence is a valid occurrence of a 0-mer, the solution now correctly outputs a dictionary. .

During development another issue that came up, was that the function did not account for all starting positions. The debugging output was helpful in identifying logical errors and ensuring that the function correctly handles various values of k and the entire sequence s.

Comparison of probability distributions¶

Now that we know how to learn probability distributions from data, we might want to compare two such distributions, for example, to test if our programs work as intended.

Let $P=\{p_1,p_2,\ldots, p_n\}$ and $Q=\{q_1,q_2,\ldots, q_n\}$ be two probability distributions for the same set of $n$ events. This means $\sum_{i=1}^n p_i=\sum_{i=1}^n q_i=1$, $0\leq p_j \leq 1$, and $0\leq q_j \leq 1$ for each event $j$.

Kullback-Leibler divergence is a measure $d()$ for the relative entropy of $P$ with respect to $Q$ defined as $d(P||Q)=\sum_{i=1}^n p_i \log\frac{p_i}{q_i}$.

This measure is always non-negative, and 0 only when $P=Q$. It can be interpreted as the gain of knowing $Q$ to encode $P$. Note that this measure is not symmetric.

19. Write function kullback_leibler to compute $d(P||Q)$. Test your solution by generating a random RNA sequence encoding the input protein sequence according to the input codon adaptation probabilities. Then you should learn the codon adaptation probabilities from the RNA sequence you generated. Then try the same with uniformly random RNA sequences (which don't have to encode any specific protein sequence). Compute the relative entropies between the three distribution (original, predicted, uniform) and you should observe a clear difference. Because $d(P||Q)$ is not symmetric, you can either print both $d(P||Q)$ and $d(Q||P)$, or their average.

This problem may be fairly tricky. Only the kullback_leibler function is automatically tested. The codon probabilities is probably a useful helper function. The main guarded section can be completed by filling out the pass sections using tooling from previous parts and fixing the placeholder lines.

Idea of solution¶

The Kullback-Leibler (KL) divergence quantifies how one probability distribution differs from a second, reference distribution. In this task, we compare the original codon distribution derived from a codon usage table, the codon distribution predicted from an RNA sequence, and a uniform distribution. The KL divergence is calculated for each pair of distributions to assess the difference in their probability distributions.

To address the task of comparing probability distributions using Kullback-Leibler (KL) divergence, we can break down the solution into several key steps:

• Generate RNA Sequences and Protein Encoding:

  • We start by generating a random protein sequence. This sequence is a string of amino acids, where each amino acid is chosen randomly from a predefined list.
  • Using the protein sequence, we then generate a corresponding RNA sequence.
  • This RNA sequence is created such that it encodes the given protein sequence according to known codon adaptation probabilities.

• Retrieve Predicted Codon Probabilities:

  • With the RNA sequence generated, we calculate the probability distribution of all possible 3-mers (codons) found within this RNA sequence. This step involves counting occurrences of each codon and normalizing these counts to obtain probabilities.
  • The function codon_probabilities() performs this calculation. It returns a dictionary where the keys are codon sequences and the values are their respective probabilities based on the RNA sequence.

• Retrieve Original Codon Probabilities:

  • We then compare our empirical probabilities with a known set of codon probabilities, obtained from a codon usage table.
  • The function get_codon_probabilities_from_table() reads the codon usage table and extracts the codon probabilities.
  • These probabilities are normalized to ensure they sum to 1.

• Retrieve Uniform Codon Probabilities:

  • We also generate a uniform probability distribution for codons.
  • In this distribution, each codon has the same probability.
  • This is done using the generate_uniform_codon_probabilities() function, which creates a uniform distribution for all possible 3-mers.

• Renormalize Distributions:

  • Before computing KL divergence, we ensure that the distributions being compared are aligned. This means filtering and renormalizing distributions to match the keys present in the reference distribution. KL divergence requires that both distributions have the same elements.
  • The renormalize() function handles this alignment, ensuring that all distributions have probabilities for the same set of codons and that these probabilities sum to 1.

• Compute Kullback-Leibler Divergence:

  • With the codon probability distributions in hand, we compute the KL divergence between different pairs of distributions: the original (reference), the predicted (empirical), and the uniform distribution (baseline).
  • The KL divergence measures how one probability distribution diverges from a second, reference probability distribution. It is calculated using the kullback_leibler() function, which sums the product of each probability in one distribution with the logarithm of the ratio of this probability to the corresponding probability in the other distribution.

By comparing these values, we can understand how well our predicted distribution matches the original codon usage and how it differs from a random distribution. Through this process, we effectively evaluate the quality of our probability predictions and assess their relevance compared to known and random baselines.

Discussion¶

The Kullback-Leibler (KL) divergence values obtained from the comparisons of the original, predicted, and uniform codon probability distributions provide insights into the quality of our predicted distribution relative to known baselines.

• A higher KL divergence value indicates a greater divergence between the two distributions, meaning more information is lost when one distribution is used to approximate the other. This suggests a poorer fit or approximation.

• A lower KL divergence value indicates that the two distributions are more similar, meaning less information is lost. This suggests a better fit or more accurate approximation.

In general, we expect the KL divergence between the original and predicted distributions to be lower than between either of them and the uniform distribution, demonstrating that specific codon usage patterns are captured better than a uniform distribution would suggest. The divergence values should clearly reflect the non-uniform nature of codon usage in real biological sequences compared to random or uniform distributions.

Here’s an analysis of the results:

• KL Divergence Between Original and Predicted Distributions:

  • d(original || predicted) = 0.15330031582470696
  • Measures how much information is lost when using the predicted distribution to approximate the original.
  • A value of 0.1533 suggests a moderate divergence, implying that while the predicted distribution is somewhat aligned with the original, there are still noticeable differences.
  • d(predicted || original) = 0.14047162716213069
  • Measures the divergence in the other direction, showing how much information is lost when the original distribution is approximated by the predicted distribution. - A value of 0.1405 is slightly lower than the previous divergence, which might indicate that the original distribution is somewhat closer to the predicted distribution than the predicted is to the original.
  • This result suggests that the predicted distribution approximates the original reasonably well but still has room for improvement.

• KL Divergence Between Original and Uniform Distributions:

  • d(original || uniform) = 0.2256062647904233
  • Shows the divergence from the original distribution to a uniform distribution, reflecting how much more structured or biased the original distribution is compared to a random distribution.
  • A value of 0.2256 suggests a substantial difference, indicating that the original codon probabilities are not uniformly distributed and are influenced by specific biological constraints or preferences.
  • d(uniform || original) = 0.09069417636570291
  • Conversely, measures how much information is lost when using a uniform distribution to approximate the original.
  • The value of 0.0907 is considerably lower than the divergence in the other direction, it indicates that the uniform distribution is not a very poor approximation of the original distribution, but it still doesn’t capture the structure of the original distribution as well as the original distribution itself.

• KL Divergence Between Predicted and Uniform Distributions:

  • (predicted || uniform) = 0.06332186437783757
  • Measures how much information is lost when approximating the uniform distribution by the predicted distribution.
  • A value of 0.0633 suggests that the predicted distribution is closer to the uniform distribution compared to the original but still retains some structure.
  • d(uniform || predicted) = 0.055725754669781385
  • Shows how much information is lost when the predicted distribution is used to approximate a uniform distribution.
  • With a value of 0.0557, this indicates that the uniform distribution is somewhat close to the predicted distribution but not as well-matched as the predicted distribution is to the uniform.

Overall, these results illustrate that while the predicted distribution captures some characteristics of the original codon usage, there is still a significant difference, indicating that the model used to generate the predicted distribution might benefit from further refinement to better approximate the complex patterns seen in the original codon probabilities.

Stationary and equilibrium distributions (extra)¶

Let us consider a Markov chain of order one on the set of nucleotides. Its transition probabilities can be expressed as a $4 \times 4$ matrix $P=(p_{ij})$, where the element $p_{ij}$ gives the probability of the $j$th nucleotide on the condition the previous nucleotide was the $i$th. An example of a transition matrix is

.

A distribution $\pi=(\pi_1,\pi_2,\pi_3,\pi_4)$ is called stationary, if $\pi = \pi P$ (the product here is matrix product).

20. Write function get_stationary_distributions that gets a transition matrix as parameter, and returns the list of stationary distributions. You can do this with NumPy by first taking transposition of both sides of the above equation to get equation $\pi^T = P^T \pi^T$. Using numpy.linalg.eig take all eigenvectors related to eigenvalue 1.0. By normalizing these vectors to sum up to one get the stationary distributions of the original transition matrix. In the main function print the stationary distributions of the above transition matrix.

Idea of solution¶

To solution aims to solve the problem of finding stationary distributions from a Markov transition matrix. Here’s a step-by-step outline:

• Transpose the Transition Matrix: To convert the problem of finding stationary distributions into one involving eigenvectors, we first transpose the transition matrix. This transformation helps us solve the problem in terms of the matrix P^T, where P is the original transition matrix.

• Compute Eigenvalues and Eigenvectors: Using NumPy's np.linalg.eig function, we compute the eigenvalues and eigenvectors of the transposed matrix. The eigenvalues provide insight into the nature of the steady-state behavior, particularly focusing on eigenvalue 1.0.

• Extract Stationary Distributions: Eigenvectors corresponding to the eigenvalue 1.0 represent the stationary distributions of the Markov chain. We normalize these eigenvectors so that their components sum to 1, converting them into valid probability distributions.

• Return the Results: The function returns these normalized eigenvectors as the list of stationary distributions.

Discussion¶

The initial function returned only one stationary distribution, likely due to filtering criteria that excluded eigenvectors with minor negative values. To address this, the filtering criteria were adjusted to account for very small negative values by treating them as zero. The function now successfully identifies stationary distributions by focusing on eigenvectors associated with the eigenvalue 1.0. This approach is appropriate for understanding the long-term behavior of the Markov chain, as stationary distributions represent the probabilities of being in each state after many transitions.

21. Implement the kl_divergence function below so that the main guarded code runs properly. Using your modified Markov chain generator generate a nucleotide sequence $s$ of length $10\;000$. Choose prefixes of $s$ of lengths $1, 10, 100, 1000$, and $10\;000$. For each of these prefixes find out their nucleotide distribution (of order 0) using your earlier tool. Use 1 as the pseudo count. Then, for each prefix, compute the KL divergence between the initial distribution and the normalized nucleotide distribution.

Idea of solution¶

NOTE. Exercises in section "Stationary and equilibrium distributions (extra)" (exercises 20, 21, and 22) are not obligatory. Thus, you only need to do 19 exercises, if you are aiming to get full points.

Discussion¶

fill in

22. Implement the following in the main function. Find the stationary distribution for the following transition matrix:

Since there is only one stationary distribution, it is called the equilibrium distribution. Choose randomly two nucleotide distributions. You can take these from your sleeve or sample them from the Dirichlet distribution. Then for each of these distributions as the initial distribution of the Markov chain, repeat the above experiment.

The main function should return tuples, where the first element is the (random) initial distribution and the second element contains the results as a list of tuples where the first element is the kl divergence and the second element the empirical nucleotide distribution, for the different prefix lengths.

The state distribution should converge to the equilibrium distribution no matter how we start the Markov chain! That is the last line of the tables should have KL-divergence very close to $0$ and an empirical distribution very close to the equilibrium distribution.

Idea of solution¶

NOTE. Exercises in section "Stationary and equilibrium distributions (extra)" (exercises 20, 21, and 22) are not obligatory. Thus, you only need to do 19 exercises, if you are aiming to get full points.

Discussion¶

fill in