<jupyter_start><jupyter_text>COMP 364 - Computer Tools for Life Sciences Homework 1Instructors: <NAME> and <NAME>**Due:  Friday, September 29th 2017 at 11:59:59 pm**In this assignment, you will learn how to implement a rudimentry profile or position weight matrix (PWM) from DNA sequences and then scan a separate DNA sequence for possible motif hits. Computational hintsPython programming visualizer: http://www.pythontutor.com/visualize.htmlmode=editAll programming concepts are drawn from lectures up to and including 'Putting it together: Nested loops and conditionals'  Biology motivation/contextFrom last week's lectures, we now know that there is a transfer of genetic information that occurs within cells from DNA-->RNA-->protein (called the central dogma of biology). In reality, this process is much more complicated as described in the figure below:If you are unfamiliar with the central dogma, we would recommend you review these documents:https://en.wikipedia.org/wiki/Central_dogma_of_molecular_biologyhttps://www.khanacademy.org/test-prep/mcat/biomolecules/amino-acids-and-proteins1/v/central-dogma-of-molecular-biology-2To help control the act of transcription (the transfer of genetic information from DNA to RNA), cells employ a set of proteins (called transcription factors) that are known to bind near genes and regulate their expression (e.g., the rate at which they are transcribed into RNA). We are able to observe the activity of transcription factors by using a technology called chromatin immunoprecipitation (or ChIP). This technology uses a targeted antibody to select for pieces of DNA that are bound by transcription factors (see figure below).The minute (or low level) details about ChIP technology are outside the scope of this course, but you must understand that the result of a ChIP protocol is a collection of DNA sequences that were bound by a specific protein of interest (i.e., a transcription factor).If you would like to learn more about ChIP technology, please see the following links:https://en.wikipedia.org/wiki/ChIP-sequencinghttp://bitesizebio.com/13541/an-introduction-to-chip-seq/Now the problem we have to address is that we have a collection of DNA sequences that were bound by a protein of interest, but we don't know the specific location where the protein attached to DNA. This is called the motif discovery problem, where:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Given a collection of DNA sequences that were bound by a protein of interest, identify the motif (or binding site) of the protein.Some of you may already be familiar with motifs and have seen examples visualized in a sequence logo:A sequence logo represents of a stack of letters (DNA nucleotides) at each position of the motif (or binding site). The relative sizes of the nucleotides indicate their frequency of representation in the collection of DNA sequences obtained from the ChIP protocol. While the total height of the letters depicts the information content of the position, in bits (a unit that you all should be familiar with now).In the upper-left panel of the figure above, we see that the binding motif of E-box is 11 nucleotides long. We see that positions 6, 7, 8, and 11 have to be nucleotides C, A, G, and G, respectively. While other positions can be multiple nucleotides or have no information about them (position 3).A profile or position weight matrix (PWM) is another way to represent the binding site of a protein that was observed to be bound to DNA. In this assignment, we will implement a PWM using only loops, lists/tuples, conditional statements, and minor arithmetic.Let's start by remembering how to declare an empty list.<jupyter_code># definition of an empty list
empty_list = [] # or list()<jupyter_output><empty_output><jupyter_text>Now, let's do something a little more useful (and biologically more relevant). Define a list object, called 'nucleotides', that contains single character representations of all possible DNA nucleotides ('A', 'C', 'G', 'T').<jupyter_code># define list of DNA nucleotides
nucleotides = [ "A", "C", "G", "T" ]<jupyter_output><empty_output><jupyter_text>With a list of nucleotides defined, we can now implement some concepts that were introduced in class. First, let's create a loop (either for or while) that iterates over each item in the 'nucleotides' list and prints the item's value & index. a) Using the built-in enumerate function ** (2.5 marks) **<jupyter_code># create a loop that iterates over each item in the list and prints out the item's index & value using enumerate()
for index, nucleotide in enumerate(nucleotides):
    print("index: ",index,"value: ",nucleotide)<jupyter_output>index:  0 value:  A
index:  1 value:  C
index:  2 value:  G
index:  3 value:  T<jupyter_text>b) Using the built-in index method ** (2.5 marks) **<jupyter_code># create a loop that iterates over each item in the list and prints out the item's index & value using index()
for nucleotide in nucleotides:
    print("index: ",nucleotides.index(nucleotide),"value: ",nucleotide)<jupyter_output>index:  0 value:  A
index:  1 value:  C
index:  2 value:  G
index:  3 value:  T<jupyter_text>Let's pretend we have just performed a ChIP protocol and observed the following six DNA sequences to be bound by our protein of interest (protein 'X'). Real ChIP-seq DNA sequences can be much longer, but for practicality, we'll limit them to jus six nucleotides.<jupyter_code># list of observed sequences where transcription factor binds
ChIP_seqs = [ "ACGTGA", "ACGGGA", "ACGCCC", "ACAGAT", "GATAGA", "ACGTTA" ]<jupyter_output><empty_output><jupyter_text>We now want to determine the frequency (how often) at which each nucleotide is bound at a given position in the motif. Remember, each sequence represents a single observation of the protein's binding motif/site. How many positions are there in the above motif? (complete the code below) ** (2.5 marks) **<jupyter_code>motif_size = len(ChIP_seqs[0])
print("motif size:",motif_size)<jupyter_output>motif size: 6<jupyter_text>For each position in the motif, we need to count the frequency of occurence for each possible nucleotide. We're going to start with just the first column and then you will implement a pair of nested loops to complete the PWM. ** (10 marks) **<jupyter_code># declare an empty list for the position weight matrix (PWM) of observed sequences
PWM = []

# define a list that contains the frequency counts for each nucleotide in the first column
nucleotide_counts = [0, 0, 0, 0] # [A, C, G, T] - note: the order is the same as our nucleotides list object above

# iterate over each DNA sequence in the list of observed DNA sequences
for DNA_seq in ChIP_seqs:
    # obtain the current observed nucleotide in the first position
    curr_nucleotide = DNA_seq[0]
    # now we can either use conditional statements or the index() function to keep track of nucleotides
    nucleotide_index = nucleotides.index(curr_nucleotide)
    # add one to the frequency count for the current nucleotide in the first position
    nucleotide_counts[nucleotide_index] += 1

# print out the nucleotide counts
print("nucleotide counts for position one:",nucleotide_counts)<jupyter_output>nucleotide counts for position one: [5, 0, 1, 0]<jupyter_text>Now we will expand our code and apply it to all positions of the motif (hint: you will need a loop). ** (10 marks) **<jupyter_code># declare an empty list for the position weight matrix (PWM) of observed sequences
PWM = []

# iterate over each position in the observed DNA sequences (hint: you can choose to use the range())
for position in range(motif_size):
    # list contains the counts of nucleotides for current column
    nucleotide_counts=[0, 0, 0, 0]
    # you can now reuse your code from the previous block (with minor adjustments)
    for DNA_seq in ChIP_seqs:
        curr_nucleotide = DNA_seq[position]
        nucleotide_index = nucleotides.index(curr_nucleotide)
        nucleotide_counts[nucleotide_index] += 1
    # append the current nucleotide counts to the PWM list object
    PWM.append(nucleotide_counts)<jupyter_output><empty_output><jupyter_text>If you have implemented the above code correctly, you should have a two-dimensional list (or a list of lists) containing frequency counts for each nucleotide at a given position in the motif. We would like to represent this in a useful table format such as:Position | A | C | G | T |:---: | :---: | :---: | :---: | :---:1 | $A_{1}$ | $C_{1}$ | $G_{1}$ | $T_{1}$ |2 | $A_{2}$ | $C_{2}$ | $G_{2}$ | $T_{2}$ |3 | $A_{3}$ | $C_{3}$ | $G_{3}$ | $T_{3}$ |4 | $A_{4}$ | $C_{4}$ | $G_{4}$ | $T_{4}$ |5 | $A_{5}$ | $C_{5}$ | $G_{5}$ | $T_{5}$ |6 | $A_{6}$ | $C_{6}$ | $G_{6}$ | $T_{6}$ |Where $A_{1}$ is the frequency count of A at position 1, $A_{2}$ is the frequency count of A at position 2, and so on.This table is called the absolute frequency matrix or position frequency matrix (PFM).To do this, we will write a loop that iterates over the two-dimensional list and prints each position's information. Note: Your table should be informative and easy to read (but it doesn't have to be perfect - i.e., don't waste your time making it look too fancy). ** (10 marks) **<jupyter_code># iterate over each position in the list object
print("\t".join(["Position","A","C","G","T"]))
for position, freqs in enumerate(PWM,1):
    print("\t".join(["\t"+str(position)]+[str(val) for val in freqs]))<jupyter_output>Position	A	C	G	T
	1	5	0	1	0
	2	1	5	0	0
	3	1	0	4	1
	4	1	1	2	2
	5	1	1	3	1
	6	4	1	0	1<jupyter_text>What should the nucleotide frequency counts in each row of your table sum to? **(2.5 marks)**6 Sequence logo creation (5 marks)Go to: http://ccg.vital-it.ch/pwmtools/pwmscan.phpFollow the instructions on the webpage and create a sequence logo of your PFM. All that should be required of you is the following:* set custom weight matrix format* paste your PFM into the text field* click submitOnce the page loads you should see a sequence logo of your PFM on the right-side of the screen. Right-click and save the image in the COMP364 HW1's subfolder called 'images'.Create a Markdown box below this text or insert beneath this text the sequence logo using the following command:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;``Replace 'your_file_name.png' with your sequence logo's filename.  Position probability matrices (PPM)Typically, PWMs are calculated as the log-likelihood ratio matrix, where $A_{1}$ represents the log-ratio of the relative frequency to the background frequency for nucleotide 'A' at position 1.First, we'll start by calculating the relative frequencies and creating the position probability matrix (PPM). To do this we will convert the absolute frequency counts to probabilities.How do we convert absolute frequencies to probabilites?Let's define some mathematical notation:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$M$ will be the PPM&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$M_{i,j}$ is the probability of nucleotide $i$ occuring at position $j$ of the PPM&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$N$ is the number of DNA sequences observed to be bound by the protein of interest&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$X_{k,j}$ is the observed nucleotide at position $j$ of observed DNA sequence $k$ By using this notation, we can define $M_{i,j}$ as a function of $X_{k,j}$ such that:$$M_{i,j}=\frac{1}{N}\sum_{k=1}^{N}I\left(X_{k,j}=i\right)$$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Where $I()$ returns 1 if $X_{k,j}=i$ and 0 otherwiseFor anyone unfamiliar with summation notation, sigma ($\sum$) is the standard for writing long sums. For example, let's sum all of the digits present in the decimal system.$$0+1+2+3+4+5+6+7+8+9=45$$This can be represented as:$$\sum_{i=0}^{9}i=0+1+2+3+4+5+6+7+8+9=45$$'$\sum_{i=0}^{9}i$' denotes that we sum all $i$ starting at $i=0$, incrememnt by 1, until we reach $i=9$.If we wanted to represent this as a for loop in Python, it would be:<jupyter_code>total = 0 # why are we not naming the variable 'sum'?
for i in range(10):
    total += i
print(total)<jupyter_output>45<jupyter_text>In practice (and as you will find), implementing the conversion from absolute frequency to a probability is relatively easy. Complete the code below to create a PPM (can be accomplished with as little as two lines of code). ** (10 marks) **Tip - use round() to limit the number of digits after the decimal point when printing the table. See the following link for more information: https://docs.python.org/3.6/library/functions.htmlround<jupyter_code># iterate over each position in the list object
N = float(motif_size)
for position, freqs in enumerate(PWM):
    PWM[position] = [freq/N for freq in freqs]

# iterate over each position in the list object
print("\t".join(["Position","A","C","G","T"]))
for position, freqs in enumerate(PWM,1):
    print("\t".join(["\t"+str(position)]+[str(round(val,3)) for val in freqs]))<jupyter_output>Position	A	C	G	T
	1	0.833	0.0	0.167	0.0
	2	0.167	0.833	0.0	0.0
	3	0.167	0.0	0.667	0.167
	4	0.167	0.167	0.333	0.333
	5	0.167	0.167	0.5	0.167
	6	0.667	0.167	0.0	0.167<jupyter_text>What should each row of your table sum to now? **(2.5 marks)**1.0 LikelihoodsNow that we have a PPM calculated, we can define the likelihood of observing a specific DNA sequences within the motif. Or in other words, if we give you a DNA sequence of 6 nucleotides, what is the probability that the protein of interest will bind it?For example, we observed the DNA sequence 'GTAGTA'. The likelihood, $L$, of of this sequence being bound by our protein of interest is:$$ L(seq\ |\ M) = M_{G,1} \times M_{T,2} \times M_{A,3} \times M_{G,4} \times M_{T,5} \times M_{A,6} $$What is the likelihood of observing the following sequences? **( 1 mark each for correct answers) ** a) CGTAGT = 0.0b) ACGTAA = 0.017146776406035666c) AAAAAA = 0.00042866941015089156(hint: fill in the code below to detemine the likelihood of observing the above sequences) ** (7 marks for correct implementation) **<jupyter_code># define sequence to calculate likelihood of
sequence = "CGTAGT"

# define initial likelihood
likelihood = 1.0

# iterate over nucleotides in sequence
for position_index, nucleotide in enumerate(sequence):
    nucleotide_index = nucleotides.index(nucleotide)
    likelihood *= PWM[position_index][nucleotide_index]
    
print("likelihood of '"+sequence+"':",likelihood)<jupyter_output>likelihood of 'CGTAGT': 0.0<jupyter_text>This is a useful program to calculate likelihoods that you have created, but it's not very practical. Add user functionality to the code: **(10 marks)*** continually prompt the user for a string containing a DNA sequence* ensure that the provided string is valid to have the likelihood calculated (print a useful error message if it does not)* then calculate the likelihood<jupyter_code>DNA_sequence_is_valid = True
while DNA_sequence_is_valid :
    # prompt user for DNA sequence
    DNA_sequence = input("Please provide a DNA sequence: ").upper().strip()
    
    # check that length of the string is equal to 6
    if( not len(DNA_sequence) == 6 ):
        print("Error - DNA sequence must be 6 nucleotides long")
        DNA_sequence_is_valid = False
    
    # check that nucleotides are valid
    for nucleotide in DNA_sequence:
        if( not nucleotide in nucleotides ):
            print("Error - Invalid nucleotide, '"+nucleotide+"' encountered")
            DNA_sequence_is_valid = False
            break
    
    if( not DNA_sequence_is_valid ):
        continue
    
    # define initial likelihood
    likelihood = 1.0

    # iterate over nucleotides in sequence
    for position_index, nucleotide in enumerate(DNA_sequence):
        nucleotide_index = nucleotides.index(nucleotide)
        likelihood *= PWM[position_index][nucleotide_index]
    
    print("likelihood of '"+DNA_sequence+"':",likelihood)<jupyter_output>Please provide a DNA sequence: ACGTAA
likelihood of 'ACGTAA': 0.017146776406035666
Please provide a DNA sequence: AAAAAA
likelihood of 'AAAAAA': 0.00042866941015089156
Please provide a DNA sequence: NNNN
Error - DNA sequence must be 6 nucleotides long
Error - Invalid nucleotide, 'N' encountered<jupyter_text>Position weight matrices (PWM)Finally, we will now convert our PPM to a PWM by calculating the log2-ratio between the observed probabilities (found in the PPM) and the background probability. For simplicity, we will assume each nucleotide occurs with an equal probability or:$$p_{i}=\frac{1}{|i|}$$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Where $|i|$ is the size of our nucleotide alphabet.What should the background frequncy be for any given nucleotide? **(2.5 marks)**$\frac{1}{4}=0.25$Now we will calculate the PWM by calculating the log2-ratio between the PPM's probabilities and the background probability, where:$$ PWM_{i,j} = log_{2}(\frac{M_{i,j}}{p_{i}})$$Again this is relatively easy and can be done in a few lines of code. **( 10 marks) **To help with the log2 transformation, we will import the Python math module (more on modules in the future) and use the math.log2() function. This function accepts a numerical object as input and returns the base-2 logarithm applied to the input. See the following link for more information: https://docs.python.org/3.4/library/math.htmlYou may encounter problems when taking the log2() of some ratios (i.e., log2(0) is negative infinity and zero divided by any number is undefined). This can be avoided by either using conditional statements or adding +1 to each ratio (feel free to choose either).<jupyter_code>import math

# iterate over each position in the list object
for position, freqs in enumerate(PWM):
    # divide by the backgound probability and calculate the log2 ratio
    PWM[position] = [math.log2((freq/0.25)+1) for freq in freqs]

# iterate over each position in the list object
print("\t".join(["Position","A","C","G","T"]))
for position, freqs in enumerate(PWM,1):
    print("\t".join(["\t"+str(position)]+[str(round(val,3)) for val in freqs]))<jupyter_output>Position	A	C	G	T
	1	2.115	0.0	0.737	0.0
	2	0.737	2.115	0.0	0.0
	3	0.737	0.0	1.874	0.737
	4	0.737	0.737	1.222	1.222
	5	0.737	0.737	1.585	0.737
	6	1.874	0.737	0.0	0.737<jupyter_text>Motif scanningAnother use for PWMs is to search for occurrences of the motif in a given DNA sequence. To do this we use what is called a sliding window approach.For example, we observe a DNA sequence of size 10 using a window of size 3.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;DNA sequence = ACGTACTATT&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Windows of size 3 = [ACG], [CGT], [GTA], [TAC], [CTA], [TAT], [ATT]In this context, the first instance of the sliding window will observe nucleotides 0-2, the second instance will then observe nucleotides 1-3, and so on. The window continues to slide along the DNA sequence until it reaches the end.For each window, we will calculate a score (similar to the likelihood) that we define as the following:$$score_{s} = \sum_{j=1}^{N}PWM_{s,j}$$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Where $s$ is the observed nucleotide in the DNA sequence located at position $j$ of the PWM. Given the defined list object below containing a random DNA sequence, identify all posible locations (above a given score threshold) where the protein may bind and sore hits based on the log2-ratio. Provide the total number of hits as well as the DNA index for the start and score of each window. ** (15 marks) **For example:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;total number of hits: ???&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;DNA index: ??	(score)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;DNA index: ??	(score)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;DNA index: ??	(score)<jupyter_code># DNA sequence to be scanned
DNA_sequence = "GTGTCGCGACAGCCCAGCATGATATCCTGAGGCGCTACGCCGACGGATGTCCCACGGAATTGCCACGGTAGCCGAGCGCTACGCGGGCGACACGAACTTATGTATGGAGCGGGCCGTCGAAAGGTCGTACCCTTGCAGTTAGCACGTAGCCCGGCCCCACTAGCACAGCAGTGCCTCGGGCGGCATCCTCATTATTAAGTTTTCTCTACAGCCAAACGACCAGGTGCGCTTCCGCGGAGCGCGGTGGAGACTCGCCCACCCGGCAGCTCTGTGACGGGGACTAAAGGGGCGATGATAATCGCGAGTGCCGCGTTATGGTGGTGTCGGGACAGAGCGGCCTTGCGGCCAGTCGTATCCCTTCTCGAGCTCCGTCCGGTTAAGCGTGACACCCCCAGCGTACCCGCAAACCGCGATGGCTGTGCTTGGGGTCGATCGCACGTAGGACGGTCCCCAACGTGAGACGCCGGGGCACCAGTTCTCACGCCCAAAGCATAAACGACGGGCAG"

# define empty list to store motif hits
hits = []
hit_threshold = 5.0

# iterate of DNA sequence and scan for possible motifs in the sequence that match our PWM
window_size = len(PWM)
for DNA_index in range(0,len(DNA_sequence)-window_size):
    window_nucleotides = DNA_sequence[DNA_index:DNA_index+window_size]
    
    # calculate the score of the window
    score = 0.0
    for position_index, nucleotide in enumerate(window_nucleotides):
        nucleotide_index = nucleotides.index(nucleotide)
        score += PWM[position_index][nucleotide_index]

    if(score > hit_threshold):
        hits.append(tuple([DNA_index,score]))
        
# print out hits to stdout
print("total number of hits:",len(hits))
for index, score in sorted(hits, key = lambda x:float(x[1]), reverse = True):
    print("DNA index: "+str(index)+"\t("+str(round(score,3))+")")<jupyter_output>total number of hits: 237
DNA index: 453	(10.787)
DNA index: 53	(9.939)
DNA index: 64	(9.939)
DNA index: 498	(9.65)
DNA index: 380	(9.409)
DNA index: 89	(9.164)
DNA index: 273	(8.913)
DNA index: 42	(8.802)
DNA index: 164	(8.802)
DNA index: 443	(8.802)
DNA index: 323	(8.672)
DNA index: 179	(8.561)
DNA index: 309	(8.561)
DNA index: 80	(8.427)
DNA index: 91	(8.316)
DNA index: 215	(8.316)
DNA index: 279	(8.316)
DNA index: 480	(8.316)
DNA index: 82	(8.271)
DNA index: 108	(8.271)
DNA index: 16	(8.271)
DNA index: 3	(8.186)
DNA index: 298	(8.186)
DNA index: 407	(8.186)
DNA index: 5	(8.075)
DNA index: 31	(8.075)
DNA index: 75	(8.075)
DNA index: 86	(8.075)
DNA index: 143	(8.065)
DNA index: 436	(8.065)
DNA index: 158	(8.027)
DNA index: 249	(8.027)
DNA index: 259	(7.824)
DNA index: 213	(7.786)
DNA index: 493	(7.786)
DNA index: 300	(7.786)
DNA index: 378	(7.786)
DNA index: 8	(7.664)
DNA index: 95	(7.664)
DNA index: 207	(7.664)
DNA index: 36	(7.579)
DNA index: 460	(7.579)
DNA index: 495	(7.579)
DNA [...]
{'fraction_non_alphanumeric': 0.0964404539446698, 'fraction_numerical': 0.03682603345335217, 'mean_word_length': 4.8356382978723405, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 55, 'https://': 5, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '606364', 'license': 'MIT', 'n_tokens_mistral': 7423, 'n_tokens_neox': 6526, 'n_words': 2863, 'path': 'hw01/A1_Solution.ipynb', 'repo_name': 'dy-lin/comp364', 'size': '21969'}