Lesson 2.6

Motif Finding Meets Oliver Cromwell

5 min read1 code challenge✓ Read

Watch the lectureHow Rolling Dice Helps Us Find Regulatory Motifs (Part 1)

What is the probability that the sun will not rise tomorrow?

In 1650, after the Scots proclaimed Charles II as king during the English Civil War, Oliver Cromwell made a famous appeal to the Church of Scotland. Urging them to see the error of their royal alliance, he pleaded,

I beseech you, in the bowels of Christ, think it possible that you may be mistaken.

The Scots rejected the appeal, and Cromwell invaded Scotland in response. His quotation later inspired the statistical maxim called Cromwell’s rule, which states that we should not use probabilities of 0 or 1 unless we are talking about logical statements that can only be true or false. In other words, we should allow a small probability for extremely unlikely events, such as "this book was written by aliens" or "the sun will not rise tomorrow". We cannot speak to the likelihood of the former event, but in the 18th Century, the French mathematician Pierre-Simon Laplace actually estimated the probability that the sun will not rise tomorrow (1/1826251), given that it has risen every day for the past 5,000 years. Although this estimate was ridiculed by his contemporaries, Laplace’s approach to this question now plays an important role in statistics.

In any observed data set, there is the possibility, especially with low-probability events or small data sets, that an event with nonzero probability does not occur. Its observed frequency is therefore zero; however, setting the empirical probability of the event equal to zero represents an inaccurate oversimplification that may cause problems. By artificially adjusting the probability of rare events, these problems can be mitigated.

Laplace’s Rule of Succession

Cromwell’s rule is relevant to the calculation of the probability of a string based on a profile matrix. For example, consider the following Profile.

The fourth symbol of TCGTGGATTTCC causes Pr(TCGTGGATTTCC, Profile) to be equal to zero. As a result, the entire string is assigned a zero probability, even though TCGTGGATTTCC differs from the consensus string at only one position. For that matter, TCGTGGATTTCC has the same low probability as AAATCTTGGAA, which is very different from the consensus string.

In order to improve this unfair scoring, bioinformaticians often substitute zeroes with small numbers called pseudocounts. The simplest approach to introducing pseudocounts, called Laplace’s Rule of Succession, is similar to the principle that Laplace used to calculate the probability that the sun will not rise tomorrow. In the case of motifs, pseudocounts often amount to adding 1 (or some other small number) to each element of Count(Motifs). For example, say we have the following motif, count, and profile matrices:

Laplace’s Rule of Succession adds 1 to each element of Count(Motifs), updating the two matrices to the following:

STOP and Think How would you use Laplace's Rule of Succession to address the shortcomings of GreedyMotifSearch?

An improved greedy motif search

The only change we need to introduce to GreedyMotifSearch in order to eliminate zeroes from the profile matrices that it constructs is to replace the sixth line of pseudocode in GreedyMotifSearch, which is highlighted below in green.

GreedyMotifSearch(Dna, k, t)
    form a set of k-mers BestMotifs by selecting 1st k-mers in each string from Dna
    for each k-mer Motif in the first string from Dna
        Motif1 ← Motif
        for i = 2 to t
            form Profile from motifs Motif1, …, Motifi - 1
            Motifi ← Profile-most probable k-mer in the i-th string in Dna
            Motifs ← (Motif1, …, Motift)
            if Score(Motifs) < Score(BestMotifs)
            BestMotifs ← Motifs
    output BestMotifs

After incorporating Laplace's Rule of Succession, GreedyMotifSearch is given by the following:

GreedyMotifSearch(Dna, k, t)
    form a set of k-mers BestMotifs by selecting 1st k-mers in each string from Dna
    for each k-mer Motif in the first string from Dna
        Motif1 ← Motif
        for i = 2 to t
            apply Laplace's Rule of Succession to form Profile from motifs Motif1, …, Motifi-1
            Motifi ← Profile-most probable k-mer in the i-th string in Dna
            Motifs ← (Motif1, …, Motift)
            if Score(Motifs) < Score(BestMotifs)
            BestMotifs ← Motifs
    output BestMotifs

We now will apply Laplace’s Rule of Succession to search for the (4, 1)-motif "ACGT" implanted in the following strings Dna:

Again, let’s assume that the algorithm has already chosen the implanted 4-mer "ACCT" from the first string. We can construct the corresponding count and profile matrices using Laplace’s Rule of Succession:

We use this profile matrix to compute the probabilities of all 4-mers in the second string from Dna:

There are two Profile-most probable 4-mers in the second string (ATGT and GTct); let’s assume that we get lucky again and choose the implanted 4-mer ATGT. We now have the following motif, count, and profile matrices:

We use this profile matrix to compute the probabilities of all 4-mers in the third string from Dna:

Again, there are two Profile-most probable 4-mers in the second string (acgG and GCGT). This time, we will assume that acgG is selected instead of the implanted 4-mer GCGT. We now have the following motif, count, and profile matrices:

We use this profile matrix to compute probabilities of all 4-mers in the fourth string from Dna:

Despite the fact that we missed the implanted 4-mer in the third string, we have now found the implanted 4-mer in the fourth string in Dna as the Profile-most probable 4-mer ACGA.

This provides us with the following motif, count, and profile matrices:

We now use this profile to compute the probabilities of all 4-mers in the fifth string in Dna:

The Profile-most probable 4-mer of the fifth string in Dna is AGGT, the implanted 4-mer. As a result, GreedyMotifSearch has produced the following motif matrix, which implies the correct consensus string ACGT:

Code Challenge Implement GreedyMotifSearch with pseudocounts. (Solve at Cogniterra or Rosalind.)Solve it with the autograderFree to start in the interactive text

• Input: Integers k and t, followed by a collection of strings Dna.
• Output: A collection of strings BestMotifs resulting from applying GreedyMotifSearch(Dna, k, t) with pseudocounts. If at any step you find more than one Profile-most probable k-mer in a given string, use the one occurring first.

You have now seen the power of pseudocounts illustrated on a small example. Running GreedyMotifSearch with pseudocounts to solve the Subtle Motif Problem returns a collection of 15-mers Motifs with Score(Motifs) = 41 and Consensus(Motifs) = AAAAAtAgaGGGGtt. Thus, Laplace's Rule of Succession has provided a great improvement over the original GreedyMotifSearch, which returned the consensus string gtAAAtAgaGatGtG with Score(Motifs) = 58.

You may be satisfied with the performance of GreedyMotifSearch, but you should know by now that we are never satisfied. Can we design an even more accurate motif finding algorithm?