Motifs

5 problems from Rosalind — Bioinformatics Textbook Track. Press Run on any block to execute it in your browser.

# Rosalind: BA2C — Find a Profile-most Probable k-mer in a String
# https://rosalind.info/problems/ba2c/
#
# Given: A string Text, an integer k, and a 4 x k matrix Profile.
# Return: A Profile-most probable k-mer in Text.

let text = "ACCTGTTTATTGCCTAAGTTCCGAACAAACCCAATATAGCCCGAGGGCCT"
let k = 5

# Rows are A, C, G, T; columns are the positions of the k-mer.
let profile = [
    [0.2, 0.2, 0.3, 0.2, 0.3],
    [0.4, 0.3, 0.1, 0.5, 0.1],
    [0.3, 0.3, 0.5, 0.2, 0.4],
    [0.1, 0.2, 0.1, 0.1, 0.2],
]

fn row_of(symbol) {
    if symbol == "A" then 0
    else if symbol == "C" then 1
    else if symbol == "G" then 2
    else 3
}

# A k-mer's probability is the product down its columns. Every entry here is
# non-zero, so no window is ruled out; a profile with a zero would silently
# eliminate one, which is why pseudocounts exist in the problems that follow.
fn probability(kmer) {
    range(0, len(kmer))
      |> reduce(|acc, i| acc * profile[row_of(substr(kmer, i, 1))][i], 1.0)
}

let kmer_windows = range(0, len(text) - k + 1) |> map(|i| substr(text, i, k))
let scores = kmer_windows |> map(|w| probability(w))
let best = kmer_windows[argmax(scores)]

println("Result:   " + best)
println("Expected: CCGAG")

fn test_ba2c_profile_most_probable() {
    assert best == "CCGAG", "BA2C: got " + best
    assert probability(best) == max(scores), "BA2C: the reported k-mer is not the most probable"
}

# Rosalind: BA2D — Implement GreedyMotifSearch
# https://rosalind.info/problems/ba2d/
#
# Given: Integers k and t, followed by a collection of strings Dna.
# Return: A collection BestMotifs resulting from GreedyMotifSearch(Dna, k, t).
# Where several profile-most probable k-mers tie, take the first.

let k = 3
let t = 5
let sequences = [
    "GGCGTTCAGGCA",
    "AAGAATCAGTCA",
    "CAAGGAGTTCGC",
    "CACGTCAATCAC",
    "CAATAATATTCG",
]

# Greedy in the strict sense: try every k-mer of the first string as a seed, then
# let each later string pick whatever its current profile likes best, never
# reconsidering. Fast, and wrong often enough that BA2E exists to fix it.
fn greedy_motifs(strings, width, pseudocount) {
    let best = strings |> map(|s| substr(s, 0, width))
    let first = strings[0]
    for start in range(0, len(first) - width + 1) {
        let motifs = [substr(first, start, width)]
        for i in range(1, len(strings)) {
            let profile = motif_profile(motifs, pseudocount)
            motifs = push(motifs, profile_most_probable(strings[i], width, profile))
        }
        if motif_score(motifs) < motif_score(best) {
            best = motifs
        }
    }
    best
}

let best_motifs = greedy_motifs(sequences, k, 0)

println("Result:")
for motif in best_motifs { println("  " + motif) }
println("Expected: CAG CAG CAA CAA CAA")

fn test_ba2d_greedy_motif_search() {
    assert join(best_motifs, " ") == "CAG CAG CAA CAA CAA",
        "BA2D: got " + join(best_motifs, " ")
    assert len(best_motifs) == t, "BA2D: one motif per string"
    # Every motif has to actually occur in its own string.
    for i in range(0, len(sequences)) {
        assert contains(sequences[i], best_motifs[i]),
            "BA2D: " + best_motifs[i] + " is not in " + sequences[i]
    }
}

# Rosalind: BA2E — Implement GreedyMotifSearch with Pseudocounts
# https://rosalind.info/problems/ba2e/
#
# Given: Integers k and t, followed by a collection of strings Dna.
# Return: BestMotifs from GreedyMotifSearch(Dna, k, t) with pseudocounts.

let k = 3
let t = 5
let sequences = [
    "GGCGTTCAGGCA",
    "AAGAATCAGTCA",
    "CAAGGAGTTCGC",
    "CACGTCAATCAC",
    "CAATAATATTCG",
]

# The same greedy loop as BA2D, changing one number. Without a pseudocount a base
# absent from a column has probability zero, so any k-mer containing it is
# impossible rather than merely unlikely — one unlucky column silently discards
# every candidate, and the search follows whichever k-mer happened to come first.
# Laplace's rule of succession fixes it by adding one to every count.
fn greedy_motifs_smoothed(strings, width, pseudocount) {
    let best = strings |> map(|s| substr(s, 0, width))
    let first = strings[0]
    for start in range(0, len(first) - width + 1) {
        let motifs = [substr(first, start, width)]
        for i in range(1, len(strings)) {
            let profile = motif_profile(motifs, pseudocount)
            motifs = push(motifs, profile_most_probable(strings[i], width, profile))
        }
        if motif_score(motifs) < motif_score(best) {
            best = motifs
        }
    }
    best
}

let best_motifs = greedy_motifs_smoothed(sequences, k, 1)

println("Result:")
for motif in best_motifs { println("  " + motif) }
println("Expected: TTC ATC TTC ATC TTC")

fn test_ba2e_greedy_motif_search_with_pseudocounts() {
    assert join(best_motifs, " ") == "TTC ATC TTC ATC TTC",
        "BA2E: got " + join(best_motifs, " ")
    for i in range(0, len(sequences)) {
        assert contains(sequences[i], best_motifs[i]),
            "BA2E: " + best_motifs[i] + " is not in " + sequences[i]
    }
    # Worth being exact about: on this five-string sample the smoothed answer
    # does *not* beat BA2D's, it ties with it — both disagree in two places. The
    # pseudocount changes which motifs are found, and pays off on inputs large
    # enough for a zero to wipe out a good candidate; a toy sample does not show
    # that, and claiming otherwise here would be checking a wish.
    let without_pseudocounts = ["CAG", "CAG", "CAA", "CAA", "CAA"]
    assert motif_score(best_motifs) == 2, "BA2E: expected a score of 2"
    assert motif_score(without_pseudocounts) == 2, "BA2D's answer also scores 2"
    assert join(best_motifs, " ") != join(without_pseudocounts, " "),
        "BA2E: the pseudocount should at least change which motifs are chosen"
}

# Rosalind: BA2F — Implement RandomizedMotifSearch
# https://rosalind.info/problems/ba2f/
#
# Given: Integers k and t, followed by a collection of strings Dna.
# Return: The best motifs found over 1000 runs of RandomizedMotifSearch.

let k = 8
let t = 5
let sequences = [
    "CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA",
    "GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG",
    "TAGTACCGAGACCGAAAGAAGTATACAGGCGT",
    "TAGATCAAGTTTCAGGTGCACGTCGGTGAACC",
    "AATCCACCAGCTCCACGTGCAATGTTGGCCTA",
]

# Seeded, so this example gives the same answer every time it is run. The
# algorithm is genuinely random; only the reporting is pinned.
set_seed(20260803)

# Start from k-mers chosen at random, then repeatedly rebuild the profile and let
# every string re-pick against it. Each round can only lower the score, so it
# stops as soon as one does not — a local optimum, and usually a poor one, which
# is why the whole thing is thrown away and restarted a thousand times.
fn one_random_run(strings, width) {
    let motifs = strings |> map(|s| {
        let start = random_int(0, len(s) - width + 1)
        substr(s, start, width)
    })
    let best = motifs
    let improving = true
    while improving {
        let profile = motif_profile(motifs, 1)
        motifs = strings |> map(|s| profile_most_probable(s, width, profile))
        if motif_score(motifs) < motif_score(best) {
            best = motifs
        } else {
            improving = false
        }
    }
    best
}

let best_motifs = one_random_run(sequences, k)
for _ in range(1, 1000) {
    let attempt = one_random_run(sequences, k)
    if motif_score(attempt) < motif_score(best_motifs) {
        best_motifs = attempt
    }
}

println("Result:")
for motif in best_motifs { println("  " + motif) }
println("Score: " + str(motif_score(best_motifs)))
println("Expected (one optimal answer): TCTCGGGG CCAAGGTG TACAGGCG TTCAGGTG TCCACGTG")
println("which also scores 9 — several distinct motif sets tie at the optimum here")

fn test_ba2f_randomized_motif_search() {
    # A randomized search is graded on the score it reaches, not on returning one
    # particular set — several sets tie at the optimum, and asserting on the
    # published one would be asserting on this seed.
    let published = ["TCTCGGGG", "CCAAGGTG", "TACAGGCG", "TTCAGGTG", "TCCACGTG"]
    assert motif_score(best_motifs) <= motif_score(published),
        "BA2F: scored " + str(motif_score(best_motifs))
            + " against the published " + str(motif_score(published))
    assert len(best_motifs) == t, "BA2F: one motif per string"
    for i in range(0, len(sequences)) {
        assert contains(sequences[i], best_motifs[i]),
            "BA2F: " + best_motifs[i] + " is not in string " + str(i)
        assert len(best_motifs[i]) == k, "BA2F: every motif is k long"
    }
}

# Rosalind: BA2G — Implement GibbsSampler
# https://rosalind.info/problems/ba2g/
#
# Given: Integers k, t and N, followed by a collection of strings Dna.
# Return: The best motifs found over 20 random starts of GibbsSampler.

let k = 8
let t = 5
let iterations = 100
let sequences = [
    "CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA",
    "GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG",
    "TAGTACCGAGACCGAAAGAAGTATACAGGCGT",
    "TAGATCAAGTTTCAGGTGCACGTCGGTGAACC",
    "AATCCACCAGCTCCACGTGCAATGTTGGCCTA",
]

set_seed(20260803)

fn all_except(items, index) {
    range(0, len(items)) |> filter(|i| i != index) |> map(|i| items[i])
}

# Not the most probable k-mer, but one drawn in proportion to its probability.
# That is the whole difference from BA2F: always taking the best makes the search
# unable to leave a local optimum, whereas sometimes taking a worse k-mer lets it
# climb back out.
fn profile_random_kmer(text, width, profile) {
    let weights = range(0, len(text) - width + 1)
        |> map(|start| profile_probability(substr(text, start, width), profile))
    let total = sum(weights)
    let target = random() * total
    let running = 0.0
    let chosen = len(weights) - 1
    for i in range(0, len(weights)) {
        running = running + weights[i]
        if running >= target and chosen == len(weights) - 1 and i < len(weights) - 1 {
            chosen = i
        }
    }
    substr(text, chosen, width)
}

# Where RandomizedMotifSearch replaces every motif at once, Gibbs replaces one at
# a time and leaves the rest standing. Changing less per step is what lets it
# keep a good partial answer instead of discarding it wholesale.
fn one_gibbs_run(strings, width, rounds) {
    let motifs = strings |> map(|s| {
        let start = random_int(0, len(s) - width + 1)
        substr(s, start, width)
    })
    let best = motifs
    for _ in range(0, rounds) {
        let i = random_int(0, len(strings))
        let profile = motif_profile(all_except(motifs, i), 1)
        motifs[i] = profile_random_kmer(strings[i], width, profile)
        if motif_score(motifs) < motif_score(best) {
            best = motifs
        }
    }
    best
}

# Rosalind suggests 20 random starts. That is a floor, not a guarantee: with 20
# this sample settles at a score of 10 rather than the optimal 9, because Gibbs
# changes one motif at a time and a bad start takes many rounds to escape. More
# starts is the knob that fixes it, and 200 finds the optimum here.
let best_motifs = one_gibbs_run(sequences, k, iterations)
for _ in range(1, 200) {
    let attempt = one_gibbs_run(sequences, k, iterations)
    if motif_score(attempt) < motif_score(best_motifs) {
        best_motifs = attempt
    }
}

println("Result:")
for motif in best_motifs { println("  " + motif) }
println("Score: " + str(motif_score(best_motifs)))
println("Expected (one optimal answer): TCTCGGGG CCAAGGTG TACAGGCG TTCAGGTG TCCACGTG, score 9")

fn test_ba2g_gibbs_sampler() {
    # Graded on the score reached, not on one particular set — as in BA2F,
    # several sets can tie. With the seed pinned above this run happens to
    # recover the published motifs exactly, which the output shows.
    let published = ["TCTCGGGG", "CCAAGGTG", "TACAGGCG", "TTCAGGTG", "TCCACGTG"]
    assert motif_score(best_motifs) <= motif_score(published),
        "BA2G: scored " + str(motif_score(best_motifs))
            + " against the published " + str(motif_score(published))
    assert len(best_motifs) == t, "BA2G: one motif per string"
    for i in range(0, len(sequences)) {
        assert contains(sequences[i], best_motifs[i]),
            "BA2G: " + best_motifs[i] + " is not in string " + str(i)
        assert len(best_motifs[i]) == k, "BA2G: every motif is k long"
    }
}