Note that RandomizedMotifSearch may change all t strings Motifs in a single iteration. This strategy may prove reckless, since some correct motifs (captured in Motifs) may potentially be discarded at the next iteration. GibbsSampler is a more cautious iterative algorithm that discards a single k-mer from the current set of motifs at each iteration and decides to either keep it or replace it with a new one. This algorithm thus moves with more caution in the space of all motifs, as illustrated below.

Randomized vs gibbs

Like RandomizedMotifSearch, GibbsSampler starts with randomly chosen k-mers in each of t DNA sequences, but it makes a random rather than a deterministic choice at each iteration. It uses randomly selected k-mers (Motif₁, …, Motif_t) to come up with another (hopefully better scoring) set of k-mers. In contrast with RandomizedMotifSearch (which deterministically defines new motifs as

Motifs(Profile(Motifs), Dna),

GibbsSampler randomly selects an integer i between 1 and t and then randomly changes only a single k-mer Motif_i.

To describe how GibbsSampler updates Motifs, we will need a slightly more advanced random number generator. Given a probability distribution (p₁, …, p_n), this random number generator, denoted Random(p₁, …, p_n), models an n-sided biased die and returns integer i with probability p_i. For example, the standard six-sided fair die represents the random number generator Random(1/6, 1/6, 1/6, 1/6, 1/6, 1/6), whereas a biased die might represent the random number generator Random(0.1, 0.2, 0.3, 0.05, 0.1, 0.25). GibbsSampler further generalizes the random number generator by using the function Random(p₁, …, p_n) defined for any set of non-negative numbers, i.e., not necessarily satisfying the condition that the p_i sum to 1. If the p_i sum to some C > 0 instead, then Random(p₁, …, p_n) is defined as Random(p₁/C, …, p_n/C), where (p₁/C, …, p_n/C) is a probability distribution. For example, for (0.1, 0.2, 0.3) with 0.1 + 0.2 + 0.3 = 0.6,

Random(0.1, 0.2, 0.3) = Random(0.1/0.6, 0.2/0.6, 0.3/0.6) = Random(1/6, 1/3, 1/2).

STOP and Think Implement the random number generator Random(p₁, …, p_n) so that it uses RandomNumber(X) (for an appropriately chosen integer X) as a subroutine.

We have previously defined the notion of a Profile-most probable k-mer in a string. We now define a Profile-randomly generated k-mer in a string Text. For each k-mer Pattern in Text, compute the probability Pr(Pattern | Profile), resulting in n = |Text| - k + 1 probabilities (p₁, …, p_n). These probabilities do not necessarily sum to 1, but we can still form the random number generator Random(p₁, …, p_n) based on them. GibbsSampler uses this random number generator to select a Profile-randomly generated k-mer at each step: if the die rolls the number i, then we define the Profile-randomly generated k-mer as the i-th k-mer in Text. While the pseudocode below repeats this procedure N times, in practice GibbsSampler depends on various stopping rules that are beyond the scope of this chapter.

GibbsSampler(Dna, k, t, N)
    randomly select k-mers Motifs = (Motif₁, …, Motif_t) in each string from Dna
    BestMotifs ← Motifs
    for j ← 1 to N
        i ← Random(t)
        Profile ← profile matrix constructed from all strings in Motifs except for Motif_i
        Motif_i ← Profile-randomly generated k-mer in the i-th sequence
        if Score(Motifs) < Score(BestMotifs)
            BestMotifs ← Motifs
    return BestMotifs

STOP and Think Note that in contrast to RandomizedMotifSearch, which always moves from higher to lower scoring Motifs, GibbsSampler may move from lower to higher scoring Motifs. Why is this reasonable?

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

Input: Integers k, t, and N, followed by a collection of strings Dna.
Output: The strings BestMotifs resulting from running GibbsSampler(Dna, k, t, N) with 20 random starts. Remember to use pseudocounts!

Note: The next lesson features a detailed example of GibbsSampler, so you may wish to return to this problem later.

Gibbs Sampling