# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# ### Counting and probability
#
# The sequence preferences of DNA-binding proteins can often be summarized as a sequence "motif" that indicates which bases are favored at each position of a binding site.
#
# One classic example of a sequence motif is the TATA-box found in many eukaryotic promoters. It's named for the first four, nearly invariant residues, but it actually extends for eight positions. The full TATA-box sequence preference can be represented as
#
# `TATA(A/T)A(A/T)(A/G)`
#
# or
#
# `TATAWAWR`
#
# using a one-letter code for nucleotide preferences.

# The TATA-box sequence motif was determined by identifying hundreds of actual TATA boxes and counting how often each nucleotide occurred.

# Here is some Python code to create a data frame with observed nucleotide counts from 389 TATA boxes taken from eukaryotic promoters (Bucher, _J Mol Biol_ (1990) 212, 563-578).
# ```
# pd.DataFrame({'A': [  16, 352,   3, 354, 268, 360, 222, 155],
#               'C': [  46,   0,  10,   0,   0,   3,   2,  44],
#               'G': [  18,   2,   2,   5,   0,  20,  44, 157],
#               'T': [ 309,  35, 374,  30, 121,   6, 121,  33]})
# ```
# Each row is a position in the TATA motif, and each column is a nucleotide.
#
# We will import the pandas package and then store this data frame in a variable called `tata_counts`.

import pandas as pd
tata_counts = pd.DataFrame({'A': [  16, 352,   3, 354, 268, 360, 222, 155],
              'C': [  46,   0,  10,   0,   0,   3,   2,  44],
              'G': [  18,   2,   2,   5,   0,  20,  44, 157],
              'T': [ 309,  35, 374,  30, 121,   6, 121,  33]})
tata_counts

# Data frames have many useful methods. For instance, we can use the .sum() method to take the sums across rows or columns. The parameter `axis=0` will calculate column sums and the parameter `axis=1` will calculate row sums. In this data frame, the row sum for each row should equal the number of sequences that Bucher examined -- each sequence will have an `A`, `C`, `G`, or `T` and so it will show up in one of the four columns.

tata_counts.sum(axis=1)

# We can then turn these counts into probabilities by dividing each nucleotide count by the total number of sequences counted. That is if 35 out of 389 TATA-box sequences have a `T` at the second position, then the probability of a `T` at position 1 in a random TATA-box sequence is 35/389, just under 10%.
#
# We'll use this to make a new `tata_probs` data frame with the _probabilities_ of each nucleotide.

tata_probs = tata_counts / 389
tata_probs

# We can now look up, e.g., the probability of a `T` at the second position, which is position 1 in Python counting
# ```
# tata_probs['T'][1]
# ```

tata_probs['T'][1]

# We're most of the way to a _probabilistic model_ of a TATA box. We will assume that each of the nucleotides in the TATA box is independent, so we can multiply these probabilities together
# $$P(\;\mathtt{TATAAAG}\;\text{for a TATA-box}\;) = 
# P(\;\mathtt{T}\mathrm{\,at\,0\;}) \times
# P(\;\mathtt{A}\mathrm{\,at\,1\;}) \times
# P(\;\mathtt{T}\mathrm{\,at\,2\;}) \times
# P(\;\mathtt{A}\mathrm{\,at\,3\;}) \times
# P(\;\mathtt{A}\mathrm{\,at\,4\;}) \times
# P(\;\mathtt{A}\mathrm{\,at\,5\;}) \times
# P(\;\mathtt{G}\mathrm{\,at\,6\;})
# $$
#
# We need to keep track of which position is which, because $P(\;\mathtt{T}\mathrm{\,at\,0\;}) \neq P(\;\mathtt{T}\mathrm{\,at\,1\;})$. The `enumerate()` function lets us keep track of a position when we're iterating over a sequence.
#
# ```
# for position, nt in enumerate(sequ):
#     print('position = ' + str(position) + ', nt = ' + str(nt))
# ```

sequ = 'TATAAAAG'
for pos, nt in enumerate(sequ):
    print('position = ' + str(pos) + ', nt = ' + str(nt))

# Now, we'll write a `for` loop to iterate over the positions in a sequence and compute a running probability.
#
# We'll start with probability 1
# ```
# prob = 1
# ```
# and then build up a running _product_ that includes the probability for each independent position
# ```
# p = tata_probs[nt][position]
# prob = prob * p
# ```
#
# We can use this to compute the probability of a "very good" TATA-box like `TATATATA`. 

prob = 1
for pos, nt in enumerate(sequ):
    prob_pos = tata_probs[nt][pos]
    prob = prob * prob_pos
    print('position ' + str(pos) + ' nt ' + str(nt) + ' prob ' + str(prob_pos))
prob

# We can also try the worst possible TATA box, `ACGCGCCT`.

sequ = 'ACGCGCCT'
prob = 1
for pos, nt in enumerate(sequ):
    prob_pos = tata_probs[nt][pos]
    prob = prob * prob_pos
    print('position ' + str(pos) + ' nt ' + str(nt) + ' prob ' + str(prob_pos))
prob

# Our final probability is 0! While $P(\;\mathtt{ACGCGCCT}\;)$ in our TATA-box model is definitely very small, it's probably not 0. We see zero `C` nucleotides at position 1 out of 389 TATA-boxes, but what if we counted 389,000? Would we find 100, 10, or 1? 
#
# We often handle these situations by adding a _pseudocount_ to our data. We add a fake count for each nucleotide, at each position, in order to eliminate zeros. The impact of this pseudocount depends on the number of real counts. If we add a pseudocount with 9 real observations, it represents 10% of our overall counts, but if we add a pseudocount with 999 real observations, it's only 0.1%.
#
# We can just add 1 to every entry to make a `tata_counts_pseudo` and use this table with pseudocounts to make our new data, and then generate a new `tata_probs` keeping in mind the "extra" counts from the pseudocounts.
#

tata_counts_pseudo = tata_counts + 1
print(tata_counts_pseudo)
print(tata_counts_pseudo.sum(axis=1))
tata_probs = tata_counts_pseudo / 393
print(tata_probs)

# Now we can use the new tata_probs to compute the probability of the best TATA-box, which is pretty similar. We can also compute the worst TATA-box, which is very low but not zero.

sequ = 'ACGCGCCT'
prob = 1
for pos, nt in enumerate(sequ):
    prob_pos = tata_probs[nt][pos]
    prob = prob * prob_pos
    print('position ' + str(pos) + ' nt ' + str(nt) + ' prob ' + str(prob_pos))
prob


# + [markdown] tags=[]
# It's getting tedious to write the same for loop every time we want to try a different sequence.
#
# We can write our own function that will compute the likelihood of a sequence under our TATA-box probability model. We'll call this function `likelihood_tata()`, for reasons that will become clearer in future lectures.
#
# ```
# def likelihood_tata(sequ):
#     prob = 1
#     for position, nt in enumerate(sequ):
#         p = tata_probs[nt][position]
#         prob = prob * p
#         print(position, nt, p, prob)
#     return(prob)
# ```
# -

def likelihood_tata(sequ):
    prob = 1
    for pos, nt in enumerate(sequ):
        p = tata_probs[nt][pos]
        prob = prob * p
    return(prob)


# Now we can easily use our function to compute the likelihood of some other possible TATA-box sequences. For example, the three sequences below are "very good" TATA-boxes that differ from the "best" TATA box at one of the three "degenerate" positions in the motif. Notice that the overall probability of getting one of these three imperfect motifs is substantially higher than the probability of the perfect TATA-box. In fact, although the TATA-box is a strong motif, fewer than 10% of the sequences generated according to our model will actually match the "best" sequence.
# ```
# TATATAAG
# TATAAATG
# TATAAAAA
# ```

likelihood_tata('TATAAAAA')

# ### _Exercise_ 
#
# We also need a probabilistic model for what _non_-TATA-box DNA sequence looks like. A model for arbitrary genomic DNA sequence is often called a "background" model.
#
# The easy background model is independent nucleotides, with probabilities determined by the overall composition of the genome. We just counted the overall number of `A`s etc in the yeast genome. A rough estimate is
#
# ```
# background = pd.Series({'A': 0.31, 'C': 0.19, 'G': 0.19, 'T': 0.31})
# ```
#
# Write a function called `likelihood_background()` to compute the probability of a sequence in the background model. Keep in mind that there is no position dependence here, and the background model is a one-dimensional Series rather than a two-dimensional DataFrame.



# Use your function to compute $P(\;\mathtt{TATAAAAG}\;)$, the probability of this very good TATA-box sequence in random genomic DNA.



# Use your function to compute $P(\;\mathtt{ACGCGCCT}\;)$, the probability of the "worst" TATA-box sequence in random genomic DNA.



StarcoderdataJupyterScriptsDedupFiltered