<jupyter_start><jupyter_text>Course Project: What's this gene?In this project you will implement a de-novo sequence assembly algorithm. You will be provided with a small sample of fragments (e.g. from an Illumina type machine) for a part of a well-known protein encoding human gene. Your task is to assemble the reads into a sequence then perform an online [BLAST](https://blast.ncbi.nlm.nih.gov/Blast.cgi) search to find out what the gene is.This project will see you use all of the techniques you will learn along the way in the course. During the course you will have the opportunity to work on the project.---  Background_Sequence Assembly_ is the process of reconstructing a sequence from small fragments (usually the output of a sequencing machine) called "**reads**". More imformation on [Wikipedia](https://en.wikipedia.org/wiki/Sequence_assembly).We will use a ["greedy algorithm"](https://en.wikipedia.org/wiki/Sequence_assemblyGreedy_algorithm) for sequence assembly. This is because it is straightforward to understand and implement but will also give you good enough results to solve the challenge. The steps of the algorithm are:1. Сalculate pairwise alignments of all fragments.1. Choose two fragments with the largest overlap.1. Merge chosen fragments.1. Repeat step 2 and 3 until only one fragment is left or you cannot merge anymore.In the following code cell, we have supplied the fragments you should assemble the sequence from.<jupyter_code>fragments = [
    "GAATTAGATAAATTCAAATTAGACTTAGGAAGGAATGTTCCCAATAGTAGACATAAAAGTCTTCGCACAGTGAAAACTAAAATGGATCAAGCAGATGATG",
    "TTCTTCAGAAGCTCCACCCTATAATTCTGAACCTGCAGAAGAATCTGAACATAAAAACAACAATTACGAACCAAACCTATTTAAAACTCCACAAAGGAAA",
    "TGCCTATTGGATCCAAAGAGAGGCCAACATTTTTTGAAATTTTTAAGACACGCTGCAACAAAGCAGATTTAGGACCAATAAGTCTTAATTGGTTTGAAGA",
    "CTCCACAAAGGAAACCATCTTATAATCAGCTGGCTTCAACTCCAATAATATTCAAAGAGCAAGGGCTGACTCTGCCGCTGTACCAATCTCCTGTAAAAGA",
    "CAACATTTTTTGAAATTTTTAAGACACGCTGCAACAAAGCAGATTTAGGACCAATAAGTCTTAATTGGTTTGAAGAACTTTCTTCAGAAGCTCCACCCTA",
    "GGCTGACTCTGCCGCTGTACCAATCTCCTGTAAAAGAATTAGATAAATTCAAATTAGACTTAGGAAGGAATGTTCCCAATAGTAGACATAAAAGTCTTCG",
    "CAAAGAGAGGCCAACATTTTTTGAAATTTTTAAGACACGCTGCAACAAAGCAGATTTAGGACCAATAAGTCTTAATTGGTTTGAAGAACTTTCTTCAGAA",
    "AGCAAGGGCTGACTCTGCCGCTGTACCAATCTCCTGTAAAAGAATTAGATAAATTCAAATTAGACTTAGGAAGGAATGTTCCCAATAGTAGACATAAAAG",
    "TCAGCTGGCTTCAACTCCAATAATATTCAAAGAGCAAGGGCTGACTCTGCCGCTGTACCAATCTCCTGTAAAAGAATTAGATAAATTCAAATTAGACTTA",
    "AATCTCCTGTAAAAGAATTAGATAAATTCAAATTAGACTTAGGAAGGAATGTTCCCAATAGTAGACATAAAAGTCTTCGCACAGTGAAAACTAAAATGGA",
    "ATGCCTATTGGATCCAAAGAGAGGCCAACATTTTTTGAAATTTTTAAGACACGCTGCAACAAAGCAGATTTAGGACCAATAAGTCTTAATTGGTTTGAAG",
    "GAATGTTCCCAATAGTAGACATAAAAGTCTTCGCACAGTGAAAACTAAAATGGATCAAGCAGATGATGTTTCCTGTCCACTTCTAAATTCTTGTCTTAGT"
]<jupyter_output><empty_output><jupyter_text>ResourcesThe `project` folder contains a complete Python solution to this problem. It also defines several _uitility_ functions that you can use when building your sequence assembly program. Please avoid looking at the solution before trying to solve each step yourself. However, when using a utility function is suggested, it will be imported and used in the template for you. --- Part 1: Getting startedYou already have enough Python knowledge to start implementing our sequence assembly program. Let's start off by breaking down the problem into more managable parts (this is a skill you will develop very quickly while programming).Step 1 says we need to compute [pairwise alignments](https://en.wikipedia.org/wiki/Sequence_alignmentPairwise_alignment). To do this we need 2 things: a way to score alignments, and a way to generate alignments. An easy way to score an alignment is known as the [_edit distance_](https://en.wikipedia.org/wiki/Levenshtein_distance), which is simply the minimum number of changes (character insertions, deletions, or substitutions) that are required to transform one string into another. Edit distance: an exampleImagine I want to find the edit distance between the strings "kitten" and "sitting":1. **k**itten → **s**itten (substitution of "s" for "k")2. sitt**e**n → sitt**i**n (substitution of "i" for "e")3. sittin → sittin**g** (insertion of "g" at the end).So the _edit distance_ between "kitten" and "sitting" is 3.Your first task is to write the beginnings of a function to compute the edit distance. Instead of a whole string, your input will be 2 characters. Please fill in the template provided for you and ensure it passes the test cases...<jupyter_code>def edit_distance(queryA, queryB):
    if _:
        return _
    
    return _

# Do not change the code below, they represent the test cases. 
assert edit_distance('A', 'T') == 1, f"I expected 1, got {edit_distance('A', 'T')}"
assert edit_distance('G', 'G') == 0, f"I expected 0, got {edit_distance('G', 'G')}"<jupyter_output><empty_output><jupyter_text>Don't worry, you will soon extend this to longer sequences but it's a good start for now.Another task we will need to do from step 3 of our algorithm is "merge" fragments. Your next task is to write a merging function that takes 2 fragments and merges them end-to-end in order. You can fill in the template below and ensure it passes the test cases...<jupyter_code>def merge(fragA, fragB):
    return _
    
assert merge("", "ATG") == "ATG", f"I expected ATG, got {merge('', 'ATG')}"
assert merge("ATG", "") == "ATG", f"I expected ATG, got {merge('ATG', '')}"
assert merge("ATG", "CCT") == "ATGCCT", f"I expected ATGCCT, got {merge('ATG', 'CCT')}"
assert merge("A", "TG") == "ATG", f"I expected ATG, got {merge('A', 'TG')}"<jupyter_output><empty_output><jupyter_text>Well done! You have now completed the beginnings of the sequence assembly program. It doesn't look like much now but it's a great foundation for tomorrow. Congratulate yourself, you've earned it!--- Part 2: Sequence assembly made easy!![How to draw an owl](https://i.kym-cdn.com/photos/images/newsfeed/000/572/078/d6d.jpg)Now that you can iterate over sequences, it's time to complete your _edit distance_ function to work on sequences longer than a single character. We will use a technique called [_dynamic programming_](https://en.wikipedia.org/wiki/Dynamic_programming). _Dynamic programming_ is a technique that allows you to solve big problems by breaking them into smaller and smaller sub-problems that are eventually "trivial". For example, the "trivial problem" when computing the edit distance is the comparison of 2 characters (that we solved yesterday). Dynamic further involves remembering (or [_memoising_](https://en.wikipedia.org/wiki/Memoization)) partial solutions as you build them up. So how is a solution "built up"? Let's have a look at the "kitten"/"sitting" example from before:Start by setting up a matrix:![Step 1](images/dp_1.jpg)Label the rows and columns:![Step 2](images/dp_2.jpg)Fill in the "boundary cases":![Step 3](images/dp_3.jpg)The general case requires you find the `min()` between three cases:![Recurrence relation](https://wikimedia.org/api/rest_v1/media/math/render/svg/10554aecc5e56da9be4657acd75b9a67b5e8b394)  More information regarding this mathematical expression is explained [here](https://en.wikipedia.org/wiki/Levenshtein_distance).If "S" is equal to "K" the substitution cost is `0`, otherwise `1`. The insertion cost and the deletion cost is always `1`, so:![General case](images/dp_4.jpg)Then systematically fill out the entire matrix:![Systematically fill in the matrix](images/dp_6.jpg)Once the matrix is complete, the bottom-left corner cell will contain the edit distance:![Edit distance](images/dp_8.jpg)I encourage you to try this for yourself with pencil and paper. Once you have convinced yourself that you understand how it works, you can implement it in Python. To begin, the memoisation matrix can be represented using a list of lists:```pythonmemo = [    -  K  I  T  T  E  N    [0, 0, 0, 0, 0, 0, 0],  -    [0, 0, 0, 0, 0, 0, 0],  S    [0, 0, 0, 0, 0, 0, 0],  I    [0, 0, 0, 0, 0, 0, 0],  T    [0, 0, 0, 0, 0, 0, 0],  T    [0, 0, 0, 0, 0, 0, 0],  I    [0, 0, 0, 0, 0, 0, 0],  N    [0, 0, 0, 0, 0, 0, 0]   G]```Lets start by working out how we read and write a value from this matrix. Let's begin by writing a function to write the value to a row and column in this matrix, and a corresponding function to get a value from a row and column in this matrix. Please fill in the functions below and ensure it passes the test cases:<jupyter_code>def set_value_in_matrix(row, column, matrix, value):
    # Set matrix at (row, column) to value
    ...
    return matrix
    
def get_value_from_matrix(row, column, matrix):
    # Get value in matrix (row, column)
    ...
    return -1
    
m = [[0, 1, 2], [3, 4, 5], [6, 7, 8]]
assert get_value_from_matrix(0, 1, set_value_in_matrix(0, 1, m, 5)) == 5
assert get_value_from_matrix(2, 0, set_value_in_matrix(2, 0, m, 81)) == 81
assert get_value_from_matrix(2, 2, set_value_in_matrix(2, 2, m, 112)) == 112<jupyter_output><empty_output><jupyter_text>Ok great work, we can read and write to our matrix. Now let's create a matrix of the correct size for our input strings. Remember that the dynamic programming matrix is the size of the input strings + 1. Please fix the function that follows and ensure it passes the tests.<jupyter_code>def create_matrix(inputA, inputB):
    row = [0] * 5
    m = []
    for _ in range(5):
        m.append(row.copy())
    
    return m

assert create_matrix('','') == [[0]], f"Expected {[[0]]}, got {create_matrix('','')}"
assert create_matrix('dn', 'a') == [[0, 0], [0, 0], [0, 0]], f"Expected {[[0, 0], [0, 0], [0, 0]]}, got {create_matrix('dn', 'a')}"<jupyter_output><empty_output><jupyter_text>Now you can initialise the matrix with the boundary cases (a edit distances against empty strings). Remember that the edit distances along the boundary are just the column or row index. For example, this is a matrix with boundary values initialised for the strings "kitten" and "sitting":```pythonmemo = [    -  K  I  T  T  E  N    [0, 1, 2, 3, 4, 5, 6],  -    [1, 0, 0, 0, 0, 0, 0],  S    [2, 0, 0, 0, 0, 0, 0],  I    [3, 0, 0, 0, 0, 0, 0],  T    [4, 0, 0, 0, 0, 0, 0],  T    [5, 0, 0, 0, 0, 0, 0],  I    [6, 0, 0, 0, 0, 0, 0],  N    [7, 0, 0, 0, 0, 0, 0]   G]```Please fix the function that follows and ensure it passes the tests.<jupyter_code>def init_matrix(inputA, inputB):
    m = create_matrix(inputA, inputB)
    rows = 0
    cols = 0
    
    for _ in range(rows):
        ...
        
    for _ in range(cols):
        ...
    
    return m

assert init_matrix('', '') == [[0]], f"Expected {[[0]]}, got {init_matrix('','')}"
assert init_matrix('dn', 'a') == [[0, 1], [1, 0], [2, 0]], f"Expected {[[0, 1], [1, 0], [2, 0]]}, got {init_matrix('dn', 'a')}"<jupyter_output><empty_output><jupyter_text>Now you're finally ready to complete the implementation of your `edit_distance()` function. You can use your `init_matrix()` function. The function template is given below, please complete and fix the function and ensure it passes the test cases.<jupyter_code>def edit_distance(queryA, queryB):
    m = init_matrix(queryA, queryB)
    rows = 0
    cols = 0

    for col in range(1, cols):
        for row in range(1, rows):
            if queryA[row - 1] == queryB[col - 1]:
                cost = 0
            else:
                cost = 1
            
            m[row][col] = min([m[row-1][col-1],
                               m[row][col - 1],
                               m[row - 1][col]
                              ])
            
    return m[0][0]

assert edit_distance('A', 'T') == 1, f"Expected 1, got {edit_distance('A', 'T')}"
assert edit_distance('G', 'G') == 0, f"Expected 0, got {edit_distance('G', 'G')}"
assert edit_distance('kitten', 'sitting') == 3, f"Expected 3, got {edit_distance('kitten', 'sitting')}"
assert edit_distance('', '') == 0, f"Expected 0, got {edit_distance('', '')}"
assert edit_distance('ABCD', 'EFGH') == 4, f"Expected 4, got {edit_distance('ABCD', 'EFGH')}"
assert edit_distance('ABCD', 'ZBCZ') == 2, f"Expected 2, got {edit_distance('ABCD', 'ZBCZ')}"<jupyter_output><empty_output><jupyter_text>Great work! You can now compute the edit distance between arbitrary strings!Now it's time to generate alignments... luckily you're already most of the way there 😊 What is an alignment?An alignment is an assignment of one-to-one correspondence between elements of a sequence. For this project you are working with sequences of nucleotides, and alignments between sequences of nucleotides are assignments of correspondence between _homologous_ nucleotides. An example of what a (pairwise) alignment looks like:            TACCCC---TCC                TA--GGTAATGGIn this example, I have highlighted "assigned correspondences" in green, insertions relative to the bottom sequence in yellow, and deletions relative to the bottom sequence in red. The `-` character is used to denote _no correspondence_. Generating alignmentsWhen you were computing the edit distance above, for each cell, you were selecting a "direction" to compute the edit distance from (up, left, or diagonal). It turns out that remembering this direction for each cell is all of the extra information you need to generate the optimal alignment for the _edit distance_ score. Then, generating the alignment is simply a matter of [_backtracking_](https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithmTracing_arrows_back_to_origin) from the bottom-right of the matrix to the top-left:![backtracking example](images/backtrack.jpg)You should try this now using pencil and paper. Create a new edit distance matrix for the strings, "saturday" and "sunday". When you fill in the matrix, draw an arrow marking the direction of the minimum value you computed for each cell. Once you finish, compare your answer with mine by clicking on the arrow below. Note that there is not a unique solution so your arrows might be slightly different to mine, but our edit distances should be the same.    Solution    Now let us see how an _alignment_ is generated using this "direction" information. Start in the bottom-right cell. If the arrow is `diagonal`, write down the column and row label for that cell. If the arrow is `left` write down the column label and a `-` (gap) character. Finally, if the arrow is `up` write down  the row label and a `-` (gap) character. Then move on the the cell pointed to by the arrow. Repeat this procedure until you arrive at the top-left cell. For example:                                                    YY                                                                            AYAY                                                                            DAYDAY                                                                            NDAYRDAY                                                                            UNDAYURDAY                                                                            -UNDAYTURDAY                                                                            --UNDAYATURDAY                                                                            S--UNDAYSATURDAY                        You should try this for yourself with pencil and paper with the strings "kitten" and "sitting".If you wish to implement this for yourself, please do as it's a great exercise to practice your Python. If you would prefer to progress without implementing this function you may use the provided `align()` function. Use of this function is demonstrated in the following code cell. Please feel free to look at the implementation and check that you understand it.<jupyter_code>from solutions.sa import align, edit_distance

assert align("sunday", "saturday", edit_distance) == ("s--unday", "saturday"), f"Expected {("s--unday", "saturday")}, got {align("sunday", "saturday", edit_distance)}"
assert align("kitten", "sitting", edit_distance) == ("kitten-", "sitting"), f"Expected {("kitten-", "sitting")}, got {align("kitten", "sitting", edit_distance)}"<jupyter_output><empty_output><jupyter_text>Scoring alignmentsAt this point you're getting very close to having all of the required building blocks to build your sequence assembly program. The final building block you require is the ability to merge fragments into longer "contigs". But we also need to be able to pick from a set of alignments, which to merge. For this we need to be able to rank alignments.For the task of sequence assembly, we expect to find _overlapping_ fragments. What does this mean? Here are some (idealised) examples of overlapping fragments:Full sequence: AAGAGGATACTATACCCTAATGAATGTTTGGAACGATTTTTCCTTCAGTGCTCTTTCAGCTCACACTCCATGGCGCATGTTGGTGACGFragment 1:        GGATACTATACCCTAATGAAFragment 2:    AAGAGGATACTATACCCTAAFragment 3:                                                                        CATGGCGCATGTTGGTGACGFragment 4:                                             CCTTCAGTGCTCTTTCAGCTFragment 5:                           ATGTTTGGAACGATTTTTCCAn alignment of `Fragment 1` and `Fragment 2` might look like this:----GGATACTATACCCTAATGAAAAGAGGATACTATACCC---T-AAThis seems like an excellent candidate to merge, but before we do that, let's look at some others. An alignment of `Fragment 1` and `Fragment 5` would look like this:-GGATACTATACCCTAATGAAATGTTTGGA-ACGATTTTTCCNot so great! What about `Fragment 2` and `Fragment 3`:-A-AGAGGATACTATACCCTAACATGGCGCAT-GT-TGGTGACGAlso not great! So how can we distinguish between a useful alignment (one that will extend 2 fragments into a contig) and an alignment that is not useful?For our purposes, _sequence identity_ will serve nicely (the supplied project code supplies a `score()` function with some slight modifications to more aggressively select useful alignments). If you wish, you may continue into the next section and use the supplied `score()` function. Otherwise, if you wish to practise here is a template for you to compute sequence identity from an alignment. Ensure your function passes the tests.<jupyter_code>def sequence_identity(top, bottom):
    identity = _
    for _ in zip(top, bottom):
        ...
        
    return identity

assert sequence_identity("", "") == 0, f"Expected 0, got {sequence_identity('', '')}"
assert sequence_identity("A", "A") == 1, f"Expected 1, got {sequence_identity('A', 'A')}"
assert sequence_identity("A", "T") == 0, f"Expected 0, got {sequence_identity('A', 'T')}"
assert sequence_identity("ATGC", "CGTA") == 0, f"Expected 0, got {sequence_identity('ATGC', 'CGTA')}"
assert sequence_identity("ATGC", "ATGC") == 4, f"Expected 4, got {sequence_identity('ATGC', 'ATGC')}"
assert sequence_identity("-AA-TG-C", "AA-TTGGG") == 3, f"Expected 3, got {sequence_identity('-AA-TG-C', 'AA-TTGGG')}"<jupyter_output><empty_output><jupyter_text>Merging alignmentsThis is the final building block you require to implement your sequence assembly program. Let's have a look at an ideal example:`----GGATACTATACCCTAATGAA``AAGAGGATACTATACCCTAA----`Excepting the gaps at the beginning and the end of the alignment, all of the nucleotides are identical. So a `merge()` for this alignment is:`AAGAGGATACTATACCCTAATGAA`Of course, it's not guaranteed that all corresponding nucleotides are identical so you will have to decide which one to accept (either top or bottom). Complete the `merge()` function below ensuring you pass the tests.<jupyter_code>def merge(top, bottom):
    ...
    
assert merge("A", "A") == "A", f"Expected 'A', got {merge('A', 'A')}"
assert merge("A", "T") == "A" or merge("A", "T") == "T", f"Expected 'A' or 'T', got {merge('A', 'T')}"
assert merge("-A", "T-") == "TA", f"Expected 'TA', got {merge('-A', 'T-')}"<jupyter_output><empty_output><jupyter_text>You have accomplished a lot today! Well done! You've implemented all of the building blocks you will require to complete the project. Tomorrow you will be able to generate sequence assemblies using what you created today, brought into existance through sheer power of your will!  Part 3: Bringing everything togetherYou have laid all the ground work. Done all of the hard work. Now all you need to do is stitch everything you've done together and produce a sequence assembly. Let's review the steps in the algorithm...1. Сalculate pairwise alignments of all fragments.1. Choose two fragments with the largest overlap.1. Merge chosen fragments.1. Repeat step 2 and 3 until only one fragment is left or you cannot merge anymore.You can translate this into a pseudocode like this:1. Loop number of fragments times    1. Calculate all pairwise alignments between fragments with `align()`    1. Pick the _best_ alignment using `score()`    1. Merge the alignment into a contig using `merge()`    1. Remove the fragments that generated the alignment and put the contig back into the collection of fragments    The supplied `project` module defines a slightly more sophisticated alignment generating algorithm that will produce better results. You're welcome to use this if you wish, it is used in the template function but if you wish to use your own alignment generation function simply remove the import.<jupyter_code>from solutions.sa import align, needleman_wunsch, score

def assemble(frags):
    n = len(frags)

    for _ in range(n - 1):
        # Generate all possible alignments
        # Generate each alignment with align(frag1, frag2, needleman_wunsch)
        
        # Pick the "best" alignment
        
        # Merge the best alignment into a contig
        contig = merge(best[0], best[1])
        
        # Remove fragments you just merged together
        
        # Add contig you just created with merge()
        
    return frags[0]

sequence = assemble(fragments.copy())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09255931812946326, 'fraction_numerical': 0.014789219073946095, 'mean_word_length': 3.8235555555555556, 'pattern_counts': {'":': 3, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 11, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 33, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '437630', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 6597, 'n_tokens_neox': 6109, 'n_words': 2604, 'path': 'projects/sequence_assembly.ipynb', 'repo_name': 'vibbits/gentle-hands-on-python', 'size': '21705'}