Motifs
6 problems from Rosalind — Bioinformatics Stronghold. Press Run on any block to execute it in your browser.
SUBS — Finding a Motif in DNA
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Finding every occurrence, overlaps included, because real motifs do overlap. The motifs worth finding are restriction sites (REVP) and binding sites — though a real binding site is a tendency rather than a fixed string, which KSIM and MPRT address.
# Rosalind: SUBS — Finding a Motif in DNA
# https://rosalind.info/problems/subs/
#
# Given: Two DNA strings s and t.
# Return: All 1-based locations of t as a substring of s, including overlaps.
let s = dna"GATATATGCATATACTT"
let t = dna"ATAT"
# find_motif reports 0-based positions; Rosalind counts from 1.
let positions = find_motif(s, t) |> map(|p| p + 1)
let result = positions |> map(|p| str(p)) |> join(" ")
println("Result: " + result)
println("Expected: 2 4 10")
fn test_subs_motif_positions() {
assert result == "2 4 10", "SUBS: got '" + result + "'"
}LCSM — Finding a Shared Motif
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Any longest common substring is a valid answer, so the assertion checks the length and that the motif really is shared, not one particular string.
# Rosalind: LCSM — Finding a Shared Motif
# https://rosalind.info/problems/lcsm/
#
# Given: A collection of DNA strings in FASTA format.
# Return: A longest common substring of the collection.
let strings = ["GATTACA", "TAGACCA", "ATACA"]
fn shared_by_all(candidate, all) {
# Bound separately: a `|> count_if(...)` directly inside a comparison gets
# read as another argument to count_if.
let hits = all |> count_if(|s| s |> contains(candidate))
hits == len(all)
}
# Search downwards from the length of the shortest string and stop at the first
# hit, so the first match found is already a longest one.
fn longest_shared(all) {
let shortest = all |> sort_by(|s| len(s)) |> first()
let width = len(shortest)
let answer = ""
while width > 0 and answer == "" {
let found = range(0, len(shortest) - width + 1)
|> map(|i| substr(shortest, i, width))
|> filter(|c| shared_by_all(c, all))
if len(found) > 0 then answer = found[0]
width = width - 1
}
answer
}
let motif = longest_shared(strings)
println("Result: " + motif + " (length " + str(len(motif)) + ")")
println("Expected: a length-2 substring shared by all, such as AC or CA or TA")
fn test_lcsm_shared_motif() {
assert len(motif) == 2, "LCSM: expected length 2, got " + str(len(motif))
assert shared_by_all(motif, strings), "LCSM: '" + motif + "' is not in every string"
}SSEQ — Finding a Spliced Motif
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A subsequence rather than a substring — the characters need not be contiguous. That is the right model for a motif split across exons, since the intervening introns are spliced out before the protein is made.
# Rosalind: SSEQ — Finding a Spliced Motif
# https://rosalind.info/problems/sseq/
#
# Given: Two DNA strings s and t.
# Return: One collection of 1-based indices of s at which t appears as a
# subsequence.
let s = "ACGTACGTGACG"
let t = "GTA"
# Greedy left-to-right: taking the earliest possible match for each character
# of t always leaves the most room for the rest.
let indices = []
let next_char = 0
let i = 0
while i < len(s) and next_char < len(t) {
if substr(s, i, 1) == substr(t, next_char, 1) {
indices = push(indices, i + 1)
next_char = next_char + 1
}
i = i + 1
}
println("Result: " + (indices |> map(|x| str(x)) |> join(" ")))
println("Expected: 3 4 5")
fn test_sseq_spliced_motif() {
assert len(indices) == len(t), "SSEQ: matched " + str(len(indices)) + " of " + str(len(t))
let spelled = indices |> map(|p| substr(s, p - 1, 1)) |> join("")
assert spelled == t, "SSEQ: indices spell '" + spelled + "'"
}MPRT — Finding a Protein Motif
solvedCLI only Problem statement Open in the workbench Download .bl
Calls NCBI, so it needs a network connection and its answer can change over time.
Fetches from UniProt, so it runs in the advisory job rather than the hermetic gate. The motif matcher itself is asserted offline — N{P}[ST]{P} has alternatives and exclusions, so it is a pattern rather than a substring search.
# Rosalind: MPRT — Finding a Protein Motif
# https://rosalind.info/problems/mprt/
#
# Given: UniProt access IDs.
# Return: For each protein containing the N-glycosylation motif, its ID and the
# 1-based positions where the motif occurs.
#
# This example fetches from UniProt, so it runs in the advisory job rather than
# the hermetic gate every other Stronghold problem passes.
let ids = ["A2Z669", "B5ZC00", "P07204_TRBM_HUMAN", "P20840_SAG1_YEAST"]
# N{P}[ST]{P} — asparagine, then anything but proline, then serine or threonine,
# then anything but proline. Not a fixed string, which is the point: a motif is a
# pattern with alternatives and exclusions, and searching for it is not
# substring matching.
fn has_motif_at(protein, i) {
substr(protein, i, 1) == "N"
and substr(protein, i + 1, 1) != "P"
and (substr(protein, i + 2, 1) == "S" or substr(protein, i + 2, 1) == "T")
and substr(protein, i + 3, 1) != "P"
}
fn motif_positions(protein) {
range(0, len(protein) - 3) |> filter(|i| has_motif_at(protein, i)) |> map(|i| i + 1)
}
# UniProt is keyed by the accession alone; the trailing name in an ID like
# P07204_TRBM_HUMAN is Rosalind's own annotation.
fn accession_of(id) { split(id, "_")[0] }
fn sequence_of(id) {
let fasta = str(uniprot_fasta(accession_of(id)))
lines(fasta) |> filter(|line| substr(line, 0, 1) != ">") |> join("")
}
let found = ids
|> map(|id| { id: id, positions: motif_positions(sequence_of(id)) })
|> filter(|entry| len(entry.positions) > 0)
println("Result:")
for entry in found {
println(" " + entry.id)
println(" " + (entry.positions |> map(|p| str(p)) |> join(" ")))
}
println("Expected: B5ZC00 85 118 142 306 395")
println(" P07204_TRBM_HUMAN 47 115 116 382 409")
println(" P20840_SAG1_YEAST 79 109 135 248 306 348 364 402 485 501 614")
fn test_mprt_finding_a_protein_motif() {
let expected = {
"B5ZC00": [85, 118, 142, 306, 395],
"P07204_TRBM_HUMAN": [47, 115, 116, 382, 409],
"P20840_SAG1_YEAST": [79, 109, 135, 248, 306, 348, 364, 402, 485, 501, 614],
}
assert len(found) == 3, "MPRT: expected 3 proteins with the motif, got " + str(len(found))
for entry in found {
assert contains(keys(expected), entry.id), "MPRT: unexpected protein " + entry.id
assert entry.positions == expected[entry.id],
"MPRT: " + entry.id + " gave " + str(entry.positions)
}
# A2Z669 has no motif, and leaving it out of the answer is part of the task.
assert (found |> count_if(|e| e.id == "A2Z669")) == 0, "MPRT: A2Z669 has no motif"
# The matcher itself, checked without the network: the exclusions are what
# make this a motif rather than a substring search.
assert motif_positions("NASA") == [1], "MPRT: NAS_ matches"
assert motif_positions("NPSA") == [], "MPRT: proline in the second position blocks it"
assert motif_positions("NASP") == [], "MPRT: proline in the fourth blocks it"
assert motif_positions("NAGA") == [], "MPRT: the third must be S or T"
assert motif_positions("NATA") == [1], "MPRT: threonine works as well as serine"
}ITWV — Finding Disjoint Motifs in a Gene
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Two motifs competing for the same characters, which is what makes this different from finding each separately. Every reachable state must consume the current character or the window ends — carrying one forward unconsumed would let the motifs match unrelated parts of the sequence.
# Rosalind: ITWV — Finding Disjoint Motifs in a Gene
# https://rosalind.info/problems/itwv/
#
# Given: A DNA string s and a collection of patterns.
# Return: The matrix M where M[j][k] = 1 if patterns j and k can be interwoven
# into some substring of s.
let text = "GACCACGGTT"
let patterns = ["ACAG", "GT", "CCG"]
# Two motifs are interwoven when a stretch of the sequence spells both at once,
# each as a subsequence, using every character for at most one of them. That is
# what makes this different from finding each motif separately: they compete for
# the same characters, so the question is whether a shared stretch can satisfy
# both — which matters when asking whether two binding sites can overlap.
#
# The state is (how much of t is placed, how much of u is placed), and the
# characters of s are consumed in order. Each character extends whichever motif
# it matches, so a cell is reachable from at most two others.
fn can_interweave(source, start, left, right) {
let rows = len(left)
let columns = len(right)
# reached[a][b] — the first a of `left` and first b of `right` are placed.
let reached = range(0, rows + 1) |> map(|_| range(0, columns + 1) |> map(|_| false))
reached[0] = reached[0]
let seed = reached[0]
seed[0] = true
reached[0] = seed
let at = start
let done = false
while at < len(source) and done == false {
let symbol = substr(source, at, 1)
let next = range(0, rows + 1) |> map(|a| range(0, columns + 1) |> map(|_| false))
for a in range(0, rows + 1) {
let line = next[a]
for b in range(0, columns + 1) {
if reached[a][b] {
# The character can be skipped only by ending the window, so
# every reachable state must consume it or the window stops
# here. Carrying it forward unconsumed is what would let the
# two motifs be found in unrelated parts of the sequence.
if a < rows and substr(left, a, 1) == symbol {
let advanced = next[a + 1]
advanced[b] = true
next[a + 1] = advanced
}
if b < columns and substr(right, b, 1) == symbol {
line[b + 1] = true
}
}
}
next[a] = line
}
reached = next
if reached[rows][columns] { done = true }
at = at + 1
}
done
}
fn interwoven_anywhere(source, left, right) {
(range(0, len(source)) |> count_if(|start| can_interweave(source, start, left, right))) > 0
}
let pairings = range(0, len(patterns)) |> map(|j|
range(0, len(patterns)) |> map(|k|
if interwoven_anywhere(text, patterns[j], patterns[k]) then 1 else 0))
println("Result:")
for line in pairings { println(" " + (line |> map(|v| str(v)) |> join(" "))) }
println("Expected: 0 0 1 / 0 1 0 / 1 0 0")
fn test_itwv_disjoint_motifs() {
let shown = pairings |> map(|line| line |> map(|v| str(v)) |> join(" ")) |> join(" / ")
assert shown == "0 0 1 / 0 1 0 / 1 0 0", "ITWV: got " + shown
# The relation is symmetric — interweaving t with u is the same question as
# interweaving u with t.
for j in range(0, len(patterns)) {
for k in range(0, len(patterns)) {
assert pairings[j][k] == pairings[k][j], "ITWV: the pairings must be symmetric"
}
}
# GT with itself works: GACCACGGTT contains GGTT, which spells GT twice using
# disjoint characters.
assert pairings[1][1] == 1, "ITWV: GT interweaves with itself"
# ACAG with itself does not — there are not enough characters to spell it
# twice disjointly anywhere.
assert pairings[0][0] == 0, "ITWV: ACAG cannot be interwoven with itself"
}KSIM — Finding All Similar Motifs
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A binding site is a tendency rather than a fixed string, so exact search finds a fraction of real sites — the motif here does not occur exactly at all. Checks every (start, length) pair directly, which is O(n^2) and stated as such: a 50 kbp genome would need the fitting-alignment form instead.
# Rosalind: KSIM — Finding All Similar Motifs
# https://rosalind.info/problems/ksim/
#
# Given: k, a motif s, and a genome t.
# Return: Every substring of t within edit distance k of s, as (position, length).
let k = 2
let motif = "ACGTAG"
let sequence = "ACGGATCGGCATCGT"
# Approximate matching, which is what real motif search always is: a binding site
# is a tendency rather than a fixed string, and an exact search finds a fraction
# of the real sites.
#
# A caveat worth stating rather than hiding: this checks every (start, length)
# pair directly, which is O(n^2) edit-distance computations. It is clear and it
# is correct, and at 15 bases it is instant — but the problem permits a 50 kbp
# genome, where it would not finish. The scalable form is a fitting alignment:
# one pass that lets the match begin anywhere in t for free, so every end
# position is scored at once instead of every pair being scored separately.
let matches = range(0, len(sequence)) |> flat_map(|start|
range(1, len(sequence) - start + 1)
|> filter(|length| edit_distance(motif, substr(sequence, start, length)) <= k)
|> map(|length| { start: start + 1, length: length }))
println("Result:")
for hit in matches { println(" " + str(hit.start) + " " + str(hit.length)) }
println("Expected: 1 4 / 1 5 / 1 6")
fn test_ksim_finding_all_similar_motifs() {
let written = matches |> map(|m| str(m.start) + " " + str(m.length))
assert sort(written) == sort(["1 4", "1 5", "1 6"]), "KSIM: got " + join(written, " / ")
# Every reported substring really is within k, and every one omitted is not —
# the second half is the part a filter can silently get wrong.
for hit in matches {
let piece = substr(sequence, hit.start - 1, hit.length)
assert edit_distance(motif, piece) <= k,
"KSIM: " + piece + " is further than " + str(k)
}
let missed = range(0, len(sequence)) |> flat_map(|start|
range(1, len(sequence) - start + 1)
|> filter(|length| edit_distance(motif, substr(sequence, start, length)) <= k)
|> filter(|length| (matches |> count_if(|m|
m.start == start + 1 and m.length == length)) == 0))
assert len(missed) == 0, "KSIM: a qualifying substring was not reported"
# An exact search finds nothing here, which is the whole point of allowing k.
assert contains(sequence, motif) == false, "KSIM: the motif does not occur exactly"
assert len(matches) > 0, "KSIM: yet three approximate matches exist"
}