Strings
29 problems from Rosalind — Bioinformatics Textbook Track. Press Run on any block to execute it in your browser.
BA1A — Compute the Number of Times a Pattern Appears in a Text
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
# Rosalind: BA1A — Compute the Number of Times a Pattern Appears in a Text
# https://rosalind.info/problems/ba1a/
#
# Given: Strings Text and Pattern.
# Return: Count(Text, Pattern), counting overlapping occurrences.
let text = "GCGCG"
let pattern = "GCG"
# find_motif reports every start position, overlaps included, so the count is
# just its length.
let result = len(find_motif(text, pattern))
println("Result: " + str(result))
println("Expected: 2")
fn test_ba1a_pattern_count() {
assert result == 2, "BA1A: got " + str(result)
}BA1B — Find the Most Frequent Words in a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Shows a real distinction: kmer_count() tallies canonical k-mers, pooling each with its reverse complement, so it reports GCAT and ATGC as one count of 4. This problem wants literal occurrences, so it counts the raw windows from kmers().
# Rosalind: BA1B — Find the Most Frequent Words in a String
# https://rosalind.info/problems/ba1b/
#
# Given: A DNA string Text and an integer k.
# Return: Every most frequent k-mer in Text.
let text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"
let k = 4
# Deliberately not kmer_count(): that tallies *canonical* k-mers, pooling each
# one with its reverse complement, which is what you want when the strand is
# unknown. Here GCAT and ATGC would merge into a single count of 4, while this
# problem asks for literal occurrences. kmers() gives the raw windows.
let window_list = kmers(dna(text), k) |> map(|km| str(km))
let distinct_count = window_list |> unique()
let tallies = distinct_count |> map(|km| { kmer: km, count: window_list |> count_if(|w| w == km) })
let best = tallies |> map(|e| e.count) |> max()
let result = tallies |> filter(|e| e.count == best) |> map(|e| e.kmer) |> sort()
println("Result: " + (result |> join(" ")) + " (each appearing " + str(best) + " times)")
println("Expected: CATG GCAT")
fn test_ba1b_most_frequent_words() {
assert len(result) == 2, "BA1B: got " + str(len(result)) + " k-mers"
assert result |> contains("CATG"), "BA1B: missing CATG"
assert result |> contains("GCAT"), "BA1B: missing GCAT"
assert best == 3, "BA1B: top count was " + str(best)
# Every window is counted exactly once.
let total = tallies |> map(|e| e.count) |> sum()
assert total == len(text) - k + 1, "BA1B: counts total " + str(total)
}BA1D — Find All Occurrences of a Pattern in a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
# Rosalind: BA1D — Find All Occurrences of a Pattern in a String
# https://rosalind.info/problems/ba1d/
#
# Given: Strings Pattern and Genome.
# Return: Every starting position where Pattern appears, zero-based.
let pattern = "ATAT"
let genome_text = "GATATATGCATATACTT"
let positions = find_motif(genome_text, pattern)
let result = positions |> map(|p| str(p)) |> join(" ")
println("Result: " + result)
println("Expected: 1 3 9")
fn test_ba1d_pattern_positions() {
assert result == "1 3 9", "BA1D: got '" + result + "'"
}BA9G — Construct the Suffix Array of a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Suffix arrays were one of three gaps this pack was built to measure. They are a builtin now, which is what LREP and MREP in the Stronghold pack stand on.
# Rosalind: BA9G — Construct the Suffix Array of a String
# https://rosalind.info/problems/ba9g/
#
# Given: A string Text.
# Return: SuffixArray(Text) — the starting positions of its suffixes, in the
# order those suffixes sort.
let text = "AACGATAGCGGTAGA$"
# The `$` is part of the input, not something the builtin adds. It sorts before
# every letter, so the empty-ish final suffix comes first, and it stops any
# suffix from being a prefix of another — which is what makes the order total.
let sa = suffix_array(text)
let result = sa |> map(|i| str(i)) |> join(", ")
println("Result: " + result)
println("Expected: 15, 14, 0, 1, 12, 6, 4, 2, 8, 13, 3, 7, 9, 10, 11, 5")
fn test_ba9g_suffix_array() {
assert result == "15, 14, 0, 1, 12, 6, 4, 2, 8, 13, 3, 7, 9, 10, 11, 5",
"BA9G: got " + result
# Every position appears exactly once, and the suffixes really are in order.
assert len(unique(sa)) == len(text), "BA9G: a position is repeated or missing"
let out_of_order = range(1, len(sa))
|> filter(|i| substr(text, sa[i - 1], len(text) - sa[i - 1])
> substr(text, sa[i], len(text) - sa[i]))
assert len(out_of_order) == 0, "BA9G: suffixes out of order at " + str(out_of_order)
}BA1E — Find Patterns Forming Clumps in a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A sliding window with a tally inside it. A k-mer qualifies once, so the answers are a set.
# Rosalind: BA1E — Find Patterns Forming Clumps in a String
# https://rosalind.info/problems/ba1e/
#
# Given: A string Genome, and integers k, L and t.
# Return: All distinct k-mers forming (L, t)-clumps in Genome.
let sequence = "CGGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAATGTACGAATGTAATGTGCCTATGGC"
let k = 5
let span = 75
let times = 4
# A k-mer forms a clump when some window of L characters contains it t times.
# Slide the window and tally what is inside it; a k-mer only needs to qualify
# once, so the answers are collected as a set.
let found = []
for start in range(0, len(sequence) - span + 1) {
let region = substr(sequence, start, span)
let counts = {}
for offset in range(0, span - k + 1) {
let mer = substr(region, offset, k)
if contains(keys(counts), mer) {
counts[mer] = counts[mer] + 1
} else {
counts[mer] = 1
}
}
for mer in keys(counts) {
if counts[mer] >= times and contains(found, mer) == false {
found = push(found, mer)
}
}
}
let result = sort(found) |> join(" ")
println("Result: " + result)
println("Expected: AATGT CGACA GAAGA (in any order)")
fn test_ba1e_clumps() {
assert len(found) == 3, "BA1E: expected 3 clumps, got " + str(len(found))
for expected in ["CGACA", "GAAGA", "AATGT"] {
assert contains(found, expected), "BA1E: missing " + expected
}
}BA1H — Find All Approximate Occurrences of a Pattern
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
hamming_distance is a builtin, so this is a windowed scan and nothing more.
# Rosalind: BA1H — Find All Approximate Occurrences of a Pattern in a String
# https://rosalind.info/problems/ba1h/
#
# Given: Strings Pattern and Text, and an integer d.
# Return: All starting positions where Pattern occurs in Text with at most d
# mismatches.
let pattern = "ATTCTGGA"
let text = "CGCCCGAATCCAGAACGCATTCCCATATTTCGGGACCACTGGCCTCCACGGTACGGACGTCAATCAAATGCCTAGCGGCTTGTGGTTTCTCCTACGCTCC"
let allowed = 3
# hamming_distance is a builtin, so the whole problem is a windowed scan: every
# position whose window is within d substitutions of the pattern.
let positions = range(0, len(text) - len(pattern) + 1)
|> filter(|i| hamming_distance(substr(text, i, len(pattern)), pattern) <= allowed)
let result = positions |> map(|i| str(i)) |> join(" ")
println("Result: " + result)
println("Expected: 6 7 26 27 78")
fn test_ba1h_approximate_matches() {
assert result == "6 7 26 27 78", "BA1H: got " + result
}BA1I — Most Frequent Words with Mismatches
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
The answer need not occur in the text at all. Tallying each window's neighbourhood is the same arithmetic over far fewer strings than testing all 4^k candidates.
# Rosalind: BA1I — Find the Most Frequent Words with Mismatches in a String
# https://rosalind.info/problems/ba1i/
#
# Given: A string Text and integers k and d.
# Return: All most frequent k-mers with up to d mismatches in Text.
let text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"
let k = 4
let d = 1
# A k-mer's count here includes every window within d substitutions of it, so
# the answer need not appear in Text at all. Rather than test every one of the
# 4^k candidates, tally the neighbourhood of each window: the same arithmetic,
# over far fewer strings.
fn neighbourhood(pattern, allowed) {
let bases = ["A", "C", "G", "T"]
let partial = [{ prefix: "", used: 0 }]
for i in range(0, len(pattern)) {
let here = substr(pattern, i, 1)
let next = []
for candidate in partial {
for base in bases {
let cost = if base == here then 0 else 1
if candidate.used + cost <= allowed {
next = push(next, { prefix: candidate.prefix + base, used: candidate.used + cost })
}
}
}
partial = next
}
partial |> map(|c| c.prefix)
}
let counts = {}
for i in range(0, len(text) - k + 1) {
for neighbour in neighbourhood(substr(text, i, k), d) {
if contains(keys(counts), neighbour) {
counts[neighbour] = counts[neighbour] + 1
} else {
counts[neighbour] = 1
}
}
}
let best = keys(counts) |> map(|mer| counts[mer]) |> max()
let winners = keys(counts) |> filter(|mer| counts[mer] == best) |> sort()
println("Result: " + (winners |> join(" ")) + " (each seen " + str(best) + " times)")
println("Expected: ATGC ATGT GATG (in any order)")
fn test_ba1i_frequent_words_with_mismatches() {
assert len(winners) == 3, "BA1I: expected 3 winners, got " + str(len(winners))
for expected in ["GATG", "ATGC", "ATGT"] {
assert contains(winners, expected), "BA1I: missing " + expected
}
}BA1J — Frequent Words with Mismatches and Reverse Complements
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A site can sit on either strand, so a k-mer is credited for its reverse complement too — which is why the winners differ from BA1I on the same input.
# Rosalind: BA1J — Find Frequent Words with Mismatches and Reverse Complements
# https://rosalind.info/problems/ba1j/
#
# Given: A DNA string Text and integers k and d.
# Return: All k-mers maximising Count_d(Text, Pattern) + Count_d(Text, revc(Pattern)).
let text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"
let k = 4
let d = 1
# Same tally as BA1I, but a k-mer is credited for its reverse complement too:
# a binding site can sit on either strand, and the counts belong together. That
# is why the winners here differ from BA1I's on the same input.
fn neighbourhood(pattern, allowed) {
let bases = ["A", "C", "G", "T"]
let partial = [{ prefix: "", used: 0 }]
for i in range(0, len(pattern)) {
let here = substr(pattern, i, 1)
let next = []
for candidate in partial {
for base in bases {
let cost = if base == here then 0 else 1
if candidate.used + cost <= allowed {
next = push(next, { prefix: candidate.prefix + base, used: candidate.used + cost })
}
}
}
partial = next
}
partial |> map(|c| c.prefix)
}
let counts = {}
for i in range(0, len(text) - k + 1) {
for neighbour in neighbourhood(substr(text, i, k), d) {
if contains(keys(counts), neighbour) {
counts[neighbour] = counts[neighbour] + 1
} else {
counts[neighbour] = 1
}
}
}
# Score each candidate together with its reverse complement.
fn tally(mer) {
if contains(keys(counts), mer) then counts[mer] else 0
}
let candidates = keys(counts)
let scores = candidates |> map(|mer| tally(mer) + tally(str(reverse_complement(dna(mer)))))
let best = max(scores)
let winners = range(0, len(candidates))
|> filter(|i| scores[i] == best)
|> map(|i| candidates[i])
|> sort()
println("Result: " + (winners |> join(" ")) + " (scoring " + str(best) + ")")
println("Expected: ACAT ATGT (in any order)")
fn test_ba1j_with_reverse_complements() {
assert len(winners) == 2, "BA1J: expected 2 winners, got " + str(len(winners))
for expected in ["ATGT", "ACAT"] {
assert contains(winners, expected), "BA1J: missing " + expected
}
# The two answers are each other's reverse complement, which is the point.
assert str(reverse_complement(dna(winners[0]))) == winners[1],
"BA1J: the winners should be a reverse-complement pair"
}BA1N — Generate the d-Neighborhood of a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Grown one position at a time, dropping any prefix that has already spent more than d substitutions.
# Rosalind: BA1N — Generate the d-Neighborhood of a String
# https://rosalind.info/problems/ba1n/
#
# Given: A DNA string Pattern and an integer d.
# Return: The collection Neighbors(Pattern, d) — every k-mer within d
# substitutions of Pattern.
let pattern = "ACG"
let d = 1
# Built by walking the pattern and, at each position, either keeping the base or
# spending one of the d substitutions on each of the other three. Written
# iteratively: grow the set of prefixes one position at a time, dropping any that
# have already overspent.
fn neighbourhood(text, allowed) {
let bases = ["A", "C", "G", "T"]
let partial = [{ prefix: "", used: 0 }]
for i in range(0, len(text)) {
let here = substr(text, i, 1)
let next = []
for candidate in partial {
for base in bases {
let cost = if base == here then 0 else 1
if candidate.used + cost <= allowed {
next = push(next, { prefix: candidate.prefix + base, used: candidate.used + cost })
}
}
}
partial = next
}
partial |> map(|c| c.prefix)
}
let result = neighbourhood(pattern, d)
println("Result: " + str(len(result)) + " neighbours")
println("Expected: 10")
println(" " + (sort(result) |> join(" ")))
fn test_ba1n_neighbourhood() {
assert len(result) == 10, "BA1N: expected 10 neighbours, got " + str(len(result))
# Every neighbour is within d, the pattern is its own neighbour, and there
# are no duplicates.
assert len(unique(result)) == len(result), "BA1N: duplicates in the neighbourhood"
assert contains(result, pattern), "BA1N: the pattern is missing from its own neighbourhood"
let too_far = result |> filter(|n| hamming_distance(n, pattern) > d)
assert len(too_far) == 0, "BA1N: these exceed d — " + str(too_far)
}BA2A — Implement MotifEnumeration
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A motif need not occur exactly anywhere, so the candidates are the neighbourhoods of the first string's windows: anything qualifying must be within d of one of them.
# Rosalind: BA2A — Implement MotifEnumeration
# https://rosalind.info/problems/ba2a/
#
# Given: Integers k and d, followed by a collection of strings Dna.
# Return: All (k, d)-motifs in Dna.
let k = 3
let d = 1
let strings = ["ATTTGGC", "TGCCTTA", "CGGTATC", "GAAAATT"]
# A (k, d)-motif appears in every string with at most d mismatches. It need not
# appear exactly anywhere, so the candidates are the neighbourhoods of the
# windows of the first string — anything qualifying must be within d of one of
# them.
fn neighbourhood(pattern, allowed) {
let bases = ["A", "C", "G", "T"]
let partial = [{ prefix: "", used: 0 }]
for i in range(0, len(pattern)) {
let here = substr(pattern, i, 1)
let next = []
for candidate in partial {
for base in bases {
let cost = if base == here then 0 else 1
if candidate.used + cost <= allowed {
next = push(next, { prefix: candidate.prefix + base, used: candidate.used + cost })
}
}
}
partial = next
}
partial |> map(|c| c.prefix)
}
fn appears_in(text, motif, allowed) {
let hits = range(0, len(text) - len(motif) + 1)
|> filter(|i| hamming_distance(substr(text, i, len(motif)), motif) <= allowed)
len(hits) > 0
}
let leader = strings[0]
let seen = []
for i in range(0, len(leader) - k + 1) {
for candidate in neighbourhood(substr(leader, i, k), d) {
if contains(seen, candidate) == false {
let missing = strings |> filter(|s| appears_in(s, candidate, d) == false)
if len(missing) == 0 {
seen = push(seen, candidate)
}
}
}
}
let result = sort(seen) |> join(" ")
println("Result: " + result)
println("Expected: ATA ATT GTT TTT")
fn test_ba2a_motif_enumeration() {
assert result == "ATA ATT GTT TTT", "BA2A: got " + result
}BA2B — Find a Median String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Every k-mer is a candidate, not only those occurring in the strings. Any minimiser is accepted, so the assertion checks the distance rather than the string: this finds ACG where the sample shows GAC, and both score 2.
# Rosalind: BA2B — Find a Median String
# https://rosalind.info/problems/ba2b/
#
# Given: An integer k and a collection of strings Dna.
# Return: A k-mer minimising d(Pattern, Dna) over all k-mers.
let k = 3
let strings = ["AAATTGACGCAT", "GACGACCACGTT", "CGTCAGCGCCTG", "GCTGAGCACCGG", "AGTACGGGACAG"]
fn distance_to(text, motif) {
range(0, len(text) - len(motif) + 1)
|> map(|i| hamming_distance(substr(text, i, len(motif)), motif))
|> min()
}
fn total_distance(motif) {
strings |> map(|s| distance_to(s, motif)) |> sum()
}
# Every k-mer is a candidate, not only those occurring in the strings — the
# median may appear in none of them. 4^k of them, enumerated as base-4 numbers.
let symbols = ["A", "C", "G", "T"]
fn number_to_pattern(index, width) {
let out = ""
let remaining = index
let position = 0
while position < width {
let power = 1
let e = 0
while e < width - position - 1 {
power = power * 4
e = e + 1
}
let digit = int(remaining / power)
out = out + symbols[digit]
remaining = remaining - digit * power
position = position + 1
}
out
}
let total_kmers = 1
let e = 0
while e < k {
total_kmers = total_kmers * 4
e = e + 1
}
let candidates = range(0, total_kmers) |> map(|i| number_to_pattern(i, k))
let scores = candidates |> map(|c| total_distance(c))
let best = min(scores)
let best_kmer = candidates[argmin(scores)]
println("Result: " + best_kmer + " (distance " + str(best) + ")")
println("Expected: GAC (any k-mer achieving the minimum is accepted)")
fn test_ba2b_median_string() {
assert total_distance(best_kmer) == best, "BA2B: the reported best_kmer is not minimal"
assert total_distance("GAC") == best, "BA2B: GAC should achieve the same minimum"
assert len(best_kmer) == k, "BA2B: wrong length"
}BA2H — Implement DistanceBetweenPatternAndStrings
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Best per string, summed across them: a motif only has to occur once in each.
# Rosalind: BA2H — Implement DistanceBetweenPatternAndStrings
# https://rosalind.info/problems/ba2h/
#
# Given: A DNA string Pattern and a collection of DNA strings Dna.
# Return: DistanceBetweenPatternAndStrings(Pattern, Dna).
let pattern = "AAA"
let strings = ["TTACCTTAAC", "GATATCTGTC", "ACGGCGTTCG", "CCCTAAAGAG", "CGTCAGAGGT"]
# The distance to one string is the best any window of it can do; the distance
# to the collection is the sum over strings. Best, not total, because a motif
# only has to occur once per string.
fn distance_to(text, motif) {
range(0, len(text) - len(motif) + 1)
|> map(|i| hamming_distance(substr(text, i, len(motif)), motif))
|> min()
}
let total = strings |> map(|s| distance_to(s, pattern)) |> sum()
println("Result: " + str(total))
println("Expected: 5")
fn test_ba2h_distance_between_pattern_and_strings() {
assert total == 5, "BA2H: got " + str(total)
# A pattern present exactly in every string would score zero.
assert distance_to("AAACCC", "AAA") == 0, "BA2H: an exact occurrence should cost nothing"
}BA9D — Find the Longest Repeat in a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
The largest entry in the LCP array, and nothing else. Two suffixes sharing a long prefix is what a repeat is, so the array has already found every one of them before the problem is read.
# Rosalind: BA9D — Find the Longest Repeat in a String
# https://rosalind.info/problems/ba9d/
#
# Given: A string Text.
# Return: A longest substring occurring more than once.
let text = "ATATCGTTTTATCGTT"
# Two suffixes sharing a long prefix is exactly a repeat, and the LCP array
# lists every such sharing between neighbours in sorted order. The longest
# repeat is the largest value in it — no search needed, only a maximum.
let padded = text + "$"
let sa = suffix_array(padded)
let lcp = lcp_array(padded)
let best = max(lcp)
let deepest = argmax(lcp)
let repeat_seq = substr(padded, sa[deepest], best)
println("Result: " + repeat_seq + " (length " + str(best) + ")")
println("Expected: TATCGTT (any longest repeat is accepted)")
fn test_ba9d_longest_repeat() {
assert len(repeat_seq) == 7, "BA9D: expected length 7, got " + str(len(repeat_seq))
# It must genuinely occur more than once.
assert len(find_motif(dna(text), dna(repeat_seq))) >= 2,
"BA9D: " + repeat_seq + " does not repeat"
assert len(find_motif(dna(text), dna("TATCGTT"))) >= 2, "BA9D: the sample answer should repeat too"
}BA9E — Find the Longest Substring Shared by Two Strings
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Concatenate with a separator, then take the largest LCP between neighbouring suffixes that came from different sides. The separator is load-bearing: without it a match can run across the join and name a substring neither string contains.
# Rosalind: BA9E — Find the Longest Substring Shared by Two Strings
# https://rosalind.info/problems/ba9e/
#
# Given: Two strings.
# Return: The longest substring occurring in both.
let first_text = "TCGGTAGATTGCGCCCACTC"
let second_text = "AGGGGCTCGCAGTGTAAGAA"
# Join the two with a separator that appears in neither, then take the suffix
# array of the pair. Neighbouring suffixes that come from *different* sides and
# share a long prefix are a shared substring — the separator is what stops a
# match running across the join and claiming something neither string contains.
let joined = first_text + "#" + second_text + "$"
let split_at = len(first_text)
let sa = suffix_array(joined)
let lcp = lcp_array(joined)
fn side_of(position) { if position < split_at then 0 else 1 }
let best = 0
let best_start = 0
for i in range(1, len(sa)) {
if side_of(sa[i]) != side_of(sa[i - 1]) and lcp[i] > best {
best = lcp[i]
best_start = sa[i]
}
}
let shared = substr(joined, best_start, best)
println("Result: " + shared + " (length " + str(best) + ")")
println("Expected: AGA (any longest shared substring is accepted)")
fn test_ba9e_longest_shared_substring() {
assert len(shared) == 3, "BA9E: expected length 3, got " + str(len(shared))
assert contains(first_text, shared), "BA9E: " + shared + " is not in the first string"
assert contains(second_text, shared), "BA9E: " + shared + " is not in the second string"
}BA9I — Construct the Burrows-Wheeler Transform of a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Sorting rotations and sorting suffixes agree once the string ends in a sentinel, so the suffix array gives the transform directly — the character before each sorted suffix.
# Rosalind: BA9I — Construct the Burrows-Wheeler Transform of a String
# https://rosalind.info/problems/ba9i/
#
# Given: A string Text.
# Return: BWT(Text).
let text = "GCGTGCCTGGTCA$"
# The BWT is the last column of the sorted rotations — but sorting rotations is
# the same as sorting suffixes when the string ends in a sentinel, so the suffix
# array gives it directly: the character just before each sorted suffix.
let n = len(text)
let sa = suffix_array(text)
let bwt = sa |> map(|i| substr(text, (i + n - 1) % n, 1)) |> join("")
println("Result: " + bwt)
println("Expected: ACTGGCT$TGCGGC")
fn test_ba9i_burrows_wheeler_transform() {
assert bwt == "ACTGGCT$TGCGGC", "BA9I: got " + bwt
# A permutation of the input, which is what makes it invertible.
assert len(bwt) == n, "BA9I: wrong length"
assert sort(chars(bwt)) == sort(chars(text)), "BA9I: not a permutation of the input"
}BA9J — Reconstruct a String from its Burrows-Wheeler Transform
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
The k-th occurrence of a symbol in the first column is the k-th in the last. That correspondence alone rebuilds the text, which is why the transform can be stored without anything beside it.
# Rosalind: BA9J — Reconstruct a String from its Burrows-Wheeler Transform
# https://rosalind.info/problems/ba9j/
#
# Given: A string Transform.
# Return: The Text whose BWT is Transform.
let transform = "TTCCTAACG$A"
# The first column is the last column sorted, and the k-th occurrence of a
# character in one column is the k-th occurrence in the other — that
# correspondence is the whole trick, and it is why the transform is reversible
# without storing anything alongside it.
let n = len(transform)
let last_column = range(0, n) |> map(|i| substr(transform, i, 1))
let first_column = sort(last_column)
# Rank each position among its own character's occurrences.
fn ranks_of(column) {
let seen = {}
let out = []
for symbol in column {
let count = if contains(keys(seen), symbol) then seen[symbol] else 0
out = push(out, count)
seen[symbol] = count + 1
}
out
}
let last_rank = ranks_of(last_column)
let first_rank = ranks_of(first_column)
# Where does row i of the last column appear in the first?
let lf = range(0, n) |> map(|i| {
let want = last_column[i]
let want_rank = last_rank[i]
(range(0, n) |> filter(|j| first_column[j] == want and first_rank[j] == want_rank))[0]
})
# Walk backwards from the row starting with the sentinel.
let walk_row = 0
let out = ""
for _ in range(0, n) {
out = last_column[walk_row] + out
walk_row = lf[walk_row]
}
# The walk produces the text rotated so the sentinel leads; move it to the end.
let text = substr(out, 1, n - 1) + substr(out, 0, 1)
println("Result: " + text)
println("Expected: TACATCACGT$")
fn test_ba9j_inverse_burrows_wheeler() {
assert text == "TACATCACGT$", "BA9J: got " + text
# Round trip: transforming the answer returns the input.
let sa = suffix_array(text)
let back = sa |> map(|i| substr(text, (i + n - 1) % n, 1)) |> join("")
assert back == transform, "BA9J: the round trip does not return the transform"
}BA9A — Construct a Trie from a Collection of Patterns
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Patterns sharing a prefix share a path, so the trie is walked once per text position no matter how many patterns there are. Node numbering is explicitly free, so the assertions are structural: every pattern spellable, no node with two edges on one symbol.
# Rosalind: BA9A — Construct a Trie from a Collection of Patterns
# https://rosalind.info/problems/ba9a/
#
# Given: A collection of strings Patterns.
# Return: The adjacency list of Trie(Patterns), each edge labelled by its symbol.
let patterns = ["ATAGA", "ATC", "GAT"]
# Patterns sharing a prefix share a path, so the trie is walked once per position
# of the text no matter how many patterns there are — which is what makes
# matching a thousand patterns cost the same as matching one. Searching each
# pattern separately would cost the text length times the pattern count.
let trie = {}
let next_id = 1
for pattern in patterns {
let node = 0
for symbol in chars(pattern) {
let existing = if contains(keys(trie), str(node)) then
trie[str(node)] |> filter(|e| e.symbol == symbol) else []
if len(existing) > 0 {
node = existing[0].to_node
} else {
let entry = { to_node: next_id, symbol: symbol }
if contains(keys(trie), str(node)) {
trie[str(node)] = push(trie[str(node)], entry)
} else {
trie[str(node)] = [entry]
}
node = next_id
next_id = next_id + 1
}
}
}
let listed = sort(keys(trie) |> map(|k| int(k))) |> flat_map(|node|
trie[str(node)] |> map(|e| str(node) + "->" + str(e.to_node) + ":" + e.symbol))
println("Result:")
for line in listed { println(" " + line) }
println("Expected: 0->1:A 1->2:T 2->3:A 3->4:G 4->5:A 2->6:C 0->7:G 7->8:A 8->9:T")
fn test_ba9a_trie_construction() {
# Node numbering is explicitly free, so the checks are structural. ATAGA and
# ATC share the prefix AT, so the trie has fewer edges than the patterns have
# characters — that saving is the whole point of a trie.
let total_characters = patterns |> map(|p| len(p)) |> sum()
assert len(listed) == 9, "BA9A: expected 9 trie, got " + str(len(listed))
assert len(listed) < total_characters,
"BA9A: a shared prefix should save trie (" + str(total_characters) + " characters)"
# Every pattern is spellable from the root, and no node has two edges on the
# same symbol.
for pattern in patterns {
let node = 0
for symbol in chars(pattern) {
let onward = trie[str(node)] |> filter(|e| e.symbol == symbol)
assert len(onward) == 1,
"BA9A: " + pattern + " has no unique path at '" + symbol + "'"
node = onward[0].to_node
}
}
for node in keys(trie) {
let symbols = trie[node] |> map(|e| e.symbol)
assert len(unique(symbols)) == len(symbols),
"BA9A: node " + node + " has two trie on one symbol"
}
}BA9B — Implement TrieMatching
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Every pattern is tested at once by a single walk, so the cost is the text length times the longest pattern rather than times the pattern count. Cross-checked against a plain scan.
# Rosalind: BA9B — Implement TrieMatching
# https://rosalind.info/problems/ba9b/
#
# Given: A string Text and a collection of strings Patterns.
# Return: Every starting position in Text where some pattern occurs.
let text = "AATCGGGTTCAATCGGGGT"
let patterns = ["ATCG", "GGGT"]
# The trie from BA9A, now used for what it is for. Starting at each position of
# the text, walk down the trie as far as the text allows; reaching a node that
# ends a pattern is a match. Every pattern is tested at once by that single walk,
# so the cost is the text length times the longest pattern, not times the number
# of patterns.
let trie = {}
let terminal = {}
let next_id = 1
for pattern in patterns {
let node = 0
for symbol in chars(pattern) {
let existing = if contains(keys(trie), str(node)) then
trie[str(node)] |> filter(|e| e.symbol == symbol) else []
if len(existing) > 0 {
node = existing[0].to_node
} else {
let entry = { to_node: next_id, symbol: symbol }
if contains(keys(trie), str(node)) {
trie[str(node)] = push(trie[str(node)], entry)
} else {
trie[str(node)] = [entry]
}
node = next_id
next_id = next_id + 1
}
}
terminal[str(node)] = true
}
let found = range(0, len(text)) |> filter(|start| {
let node = 0
let at = start
let matched = false
let walking = true
while walking {
if contains(keys(terminal), str(node)) {
matched = true
walking = false
} else {
if at >= len(text) or contains(keys(trie), str(node)) == false {
walking = false
} else {
let onward = trie[str(node)] |> filter(|e| e.symbol == substr(text, at, 1))
if len(onward) == 0 {
walking = false
} else {
node = onward[0].to_node
at = at + 1
}
}
}
}
matched
})
println("Result: " + (found |> map(|p| str(p)) |> join(" ")))
println("Expected: 1 4 11 15")
fn test_ba9b_trie_matching() {
assert (found |> map(|p| str(p)) |> join(" ")) == "1 4 11 15",
"BA9B: got " + str(found)
# Every reported position really starts one of the patterns...
for position in found {
let here = patterns |> filter(|p| substr(text, position, len(p)) == p)
assert len(here) > 0, "BA9B: nothing matches at " + str(position)
}
# ...and nothing was missed, checked the slow way against the trie's answer.
let by_brute_force = range(0, len(text))
|> filter(|start| (patterns |> count_if(|p| substr(text, start, len(p)) == p)) > 0)
assert found == by_brute_force,
"BA9B: the trie found " + str(found) + " but scanning finds " + str(by_brute_force)
}BA9H — Pattern Matching with the Suffix Array
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
The same answer as BA9B reached the other way, and the example says which to reach for: a trie is built from the patterns and suits many patterns against one text; a suffix array is built from the text and suits one text queried repeatedly.
# Rosalind: BA9H — Pattern Matching with the Suffix Array
# https://rosalind.info/problems/ba9h/
#
# Given: A string Text and a collection of strings Patterns.
# Return: Every starting position in Text where some pattern occurs.
let text = "AATCGGGTTCAATCGGGGT"
let patterns = ["ATCG", "GGGT"]
# The same answer as BA9B, reached the other way. Every occurrence of a pattern
# is the start of a suffix beginning with it, and the suffix array holds the
# suffixes in sorted order — so those suffixes form one contiguous band, findable
# by binary search.
#
# The trade against the trie: the trie is built from the patterns and scans the
# text once, so it suits many patterns against one text. The suffix array is
# built from the text and searched per pattern, so it suits one text queried
# repeatedly — a reference genome, indexed once.
let sa = suffix_array(text)
fn suffix_at(source, start) { substr(source, start, len(source) - start) }
fn starts_with(haystack, needle) {
len(haystack) >= len(needle) and substr(haystack, 0, len(needle)) == needle
}
fn band_for(pattern, order, source) {
# Lower edge: the first suffix not less than the pattern.
let low = 0
let high = len(order)
while low < high {
let middle = floor((low + high) / 2)
if suffix_at(source, order[middle]) < pattern { low = middle + 1 } else { high = middle }
}
let start = low
# Upper edge: the first suffix that no longer begins with it.
high = len(order)
while low < high {
let middle = floor((low + high) / 2)
if starts_with(suffix_at(source, order[middle]), pattern) {
low = middle + 1
} else {
high = middle
}
}
{ start: start, stop: low }
}
let found = sort(patterns |> flat_map(|pattern| {
let band = band_for(pattern, sa, text)
range(band.start, band.stop) |> map(|i| sa[i])
}))
println("Result: " + (found |> map(|p| str(p)) |> join(" ")))
println("Expected: 1 4 11 15")
fn test_ba9h_suffix_array_matching() {
assert (found |> map(|p| str(p)) |> join(" ")) == "1 4 11 15", "BA9H: got " + str(found)
# Same answer as scanning, which is the only thing that matters.
let by_brute_force = sort(patterns |> flat_map(|p|
range(0, len(text) - len(p) + 1) |> filter(|s| substr(text, s, len(p)) == p)))
assert found == by_brute_force,
"BA9H: the suffix array found " + str(found) + " but scanning finds " + str(by_brute_force)
# A pattern that is not there returns an empty band rather than a wrong one.
let missing = band_for("TTTT", sa, text)
assert missing.stop == missing.start, "BA9H: TTTT should match nothing"
}BA9K — Generate the Last-to-First Mapping of a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
The kth occurrence of a symbol in one column is the kth in the other, so a row is found by counting rather than searching. Asserted to be a bijection over every row, which is what makes the walk in BA9J terminate.
# Rosalind: BA9K — Generate the Last-to-First Mapping of a String
# https://rosalind.info/problems/ba9k/
#
# Given: A string Transform and an integer i.
# Return: LastToFirst(i) — where the symbol at position i of the last column sits
# in the first column.
let transform = "T$GACCA"
let position = 3
# The first column is the last column sorted, and — this is the property
# everything else rests on — the kth occurrence of a symbol in one column is the
# kth occurrence in the other. So a symbol's row in the first column is found by
# counting how many of its own kind precede it in the last, not by searching.
#
# BA9J used this to undo the transform; BA9L uses it to search without ever
# rebuilding the text.
let last_column = chars(transform)
let first_column = sort(last_column)
fn occurrence_rank(column, at) {
range(0, at) |> count_if(|i| column[i] == column[at])
}
let symbol = last_column[position]
let which_occurrence = occurrence_rank(last_column, position)
# That same occurrence of the symbol, located in the first column.
let matching = range(0, len(first_column)) |> filter(|i| first_column[i] == symbol)
let answer = matching[which_occurrence]
println("Result: " + str(answer))
println("Expected: 1")
fn test_ba9k_last_to_first() {
assert answer == 1, "BA9K: got " + str(answer)
assert first_column[answer] == symbol, "BA9K: the two columns must agree on the symbol"
# The mapping is a bijection: every row of the last column lands on a
# different row of the first. That is what makes the walk in BA9J terminate.
let images = range(0, len(last_column)) |> map(|i| {
let s = last_column[i]
let r = occurrence_rank(last_column, i)
(range(0, len(first_column)) |> filter(|j| first_column[j] == s))[r]
})
assert len(unique(images)) == len(last_column), "BA9K: last-to-first must be a bijection"
assert sort(images) == range(0, len(last_column)), "BA9K: and must cover every row"
}BA9L — Implement BWMatching
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Searching the transform without ever rebuilding the text: sorted rows mean every match forms one contiguous band, narrowed one symbol at a time. Cross-checked by inverting the transform — the point being that the search never needed to.
# Rosalind: BA9L — Implement BWMatching
# https://rosalind.info/problems/ba9l/
#
# Given: A string BWT(Text) and a collection of Patterns.
# Return: How many times each pattern occurs in Text.
let transform = "TCCTCTATGAGATCCTATTCTATGAAACCTTCA$GACCAAAATTCTCCGGC"
let patterns = ["CCT", "CAC", "GAG", "CAG", "ATC"]
# Searching the transform directly, without ever rebuilding the text. Rows of the
# Burrows-Wheeler matrix are sorted, so every row beginning with a given suffix
# forms one contiguous band. Reading the pattern backwards, each new symbol
# narrows the band by one last-to-first step, and the band's width at the end is
# the number of occurrences.
#
# This is why a read aligner can index a genome once and answer queries against
# the index alone: the cost depends on the pattern, not on the genome.
let last_column = chars(transform)
let first_column = sort(last_column)
fn occurrence_rank(column, at) {
range(0, at) |> count_if(|i| column[i] == column[at])
}
let last_to_first = range(0, len(last_column)) |> map(|i| {
let symbol = last_column[i]
let rank = occurrence_rank(last_column, i)
(range(0, len(first_column)) |> filter(|j| first_column[j] == symbol))[rank]
})
fn count_matches(pattern, last, mapping) {
let top = 0
let bottom = len(last) - 1
let remaining = reverse(pattern)
let searching = true
for symbol in chars(remaining) {
if searching {
let rows = range(top, bottom + 1) |> filter(|i| last[i] == symbol)
if len(rows) == 0 {
top = 1
bottom = 0
searching = false
} else {
top = mapping[rows[0]]
bottom = mapping[rows[len(rows) - 1]]
}
}
}
if bottom >= top then bottom - top + 1 else 0
}
let counts = patterns |> map(|p| count_matches(p, last_column, last_to_first))
println("Result: " + (counts |> map(|c| str(c)) |> join(" ")))
println("Expected: 2 1 1 0 1")
fn test_ba9l_bw_matching() {
assert (counts |> map(|c| str(c)) |> join(" ")) == "2 1 1 0 1",
"BA9L: got " + str(counts)
# Checked against the text itself, recovered by inverting the transform as in
# BA9J — the point being that BWMatching never needed it.
let n = len(last_column)
let out = ""
let walk_row = 0
for _ in range(0, n) {
out = last_column[walk_row] + out
walk_row = last_to_first[walk_row]
}
let text = substr(out, 1, n - 1) + substr(out, 0, 1)
for i in range(0, len(patterns)) {
let by_scanning = range(0, len(text) - len(patterns[i]) + 1)
|> count_if(|s| substr(text, s, len(patterns[i])) == patterns[i])
assert counts[i] == by_scanning,
"BA9L: " + patterns[i] + " counted " + str(counts[i])
+ " but occurs " + str(by_scanning) + " times"
}
assert counts[3] == 0, "BA9L: CAG does not occur"
}BA9M — Implement BetterBWMatching
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Two precomputed tables remove BWMatching's scan, so a query costs time proportional to the pattern rather than to the text. That difference is what makes indexing a genome once and querying it billions of times practical.
# Rosalind: BA9M — Implement BetterBWMatching
# https://rosalind.info/problems/ba9m/
#
# Given: A string BWT(Text) and a collection of Patterns.
# Return: How many times each pattern occurs in Text.
let transform = "GGCGCCGC$TAGTCACACACGCCGTA"
let patterns = ["ACC", "CCG", "CAG"]
# The same search as BA9L, made fast. BWMatching scans the band on every step to
# find the rows carrying the next symbol, so a query costs time proportional to
# the text. Precomputing two small tables removes that scan entirely:
#
# first_occurrence where each symbol's block begins in the first column
# occurrences_before how many of each symbol appear in the last column before
# each row
#
# With those, narrowing the band is two lookups and an addition, so a query costs
# time proportional to the *pattern*. That difference is what makes indexing a
# genome once and querying it billions of times practical.
let last_column = chars(transform)
let first_column = sort(last_column)
let alphabet = sort(unique(last_column))
let first_occurrence = {}
for symbol in alphabet {
first_occurrence[symbol] = (range(0, len(first_column))
|> filter(|i| first_column[i] == symbol))[0]
}
# count[symbol][i] = occurrences of symbol in the first i entries of the last
# column.
let occurrences_before = {}
for symbol in alphabet {
let running = [0]
for i in range(0, len(last_column)) {
running = push(running, running[i] + (if last_column[i] == symbol then 1 else 0))
}
occurrences_before[symbol] = running
}
fn count_matches(pattern, alpha, firsts, counts, length) {
let top = 0
let bottom = length - 1
let searching = true
for symbol in chars(reverse(pattern)) {
if searching {
if contains(alpha, symbol) == false {
searching = false
top = 1
bottom = 0
} else {
let before = counts[symbol][top]
let within = counts[symbol][bottom + 1]
if within - before == 0 {
searching = false
top = 1
bottom = 0
} else {
top = firsts[symbol] + before
bottom = firsts[symbol] + within - 1
}
}
}
}
if bottom >= top then bottom - top + 1 else 0
}
let counts = patterns |> map(|p|
count_matches(p, alphabet, first_occurrence, occurrences_before, len(last_column)))
println("Result: " + (counts |> map(|c| str(c)) |> join(" ")))
println("Expected: 1 2 1")
fn test_ba9m_better_bw_matching() {
assert (counts |> map(|c| str(c)) |> join(" ")) == "1 2 1", "BA9M: got " + str(counts)
# Checked against the text, recovered by inverting the transform.
let n = len(last_column)
let mapping = range(0, n) |> map(|i| {
let symbol = last_column[i]
let rank = range(0, i) |> count_if(|j| last_column[j] == symbol)
first_occurrence[symbol] + rank
})
let out = ""
let walk_row = 0
for _ in range(0, n) {
out = last_column[walk_row] + out
walk_row = mapping[walk_row]
}
let text = substr(out, 1, n - 1) + substr(out, 0, 1)
for i in range(0, len(patterns)) {
let by_scanning = range(0, len(text) - len(patterns[i]) + 1)
|> count_if(|s| substr(text, s, len(patterns[i])) == patterns[i])
assert counts[i] == by_scanning,
"BA9M: " + patterns[i] + " counted " + str(counts[i])
+ " but occurs " + str(by_scanning) + " times"
}
}BA9Q — Construct the Partial Suffix Array of a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A full suffix array of a human genome is 12 GB before the sequence itself. Keeping every Kth value cuts that by a factor of K, with the rest recoverable by walking the BWT — the compromise real read aligners ship with.
# Rosalind: BA9Q — Construct the Partial Suffix Array of a String
# https://rosalind.info/problems/ba9q/
#
# Given: A string Text and a positive integer K.
# Return: The partial suffix array — the pairs (i, SuffixArray(i)) whose value is
# a multiple of K.
let text = "PANAMABANANAS$"
let step = 5
# A full suffix array stores one integer per position, which for a human genome
# is 12 GB before the sequence itself. Keeping only every Kth *value* cuts that by
# a factor of K, and the discarded entries are recoverable by walking backwards
# through the BWT until a kept one is reached — trading a little time per query
# for memory that would otherwise not fit. This is the compromise real
# read aligners like Bowtie ship with.
let full = suffix_array(text)
let partial = range(0, len(full))
|> filter(|i| full[i] % step == 0)
|> map(|i| { row: i, value: full[i] })
println("Result:")
for entry in partial { println(" " + str(entry.row) + "," + str(entry.value)) }
println("Expected: 1,5 / 11,10 / 12,0")
fn test_ba9q_partial_suffix_array() {
let written = partial |> map(|e| str(e.row) + "," + str(e.value))
assert join(written, " / ") == "1,5 / 11,10 / 12,0", "BA9Q: got " + join(written, " / ")
# Every kept value is a multiple of K, and every multiple present in the full
# array is kept — nothing is dropped that should not be.
for entry in partial {
assert entry.value % step == 0, "BA9Q: " + str(entry.value) + " is not a multiple of 5"
assert full[entry.row] == entry.value, "BA9Q: row " + str(entry.row) + " disagrees"
}
let expected_count = full |> count_if(|v| v % step == 0)
assert len(partial) == expected_count, "BA9Q: some multiples were dropped"
# The saving: three entries kept out of fourteen.
assert len(partial) < len(full), "BA9Q: a partial array must be smaller"
}BA9C — Construct the Suffix Tree of a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
A trie of every suffix with non-branching chains collapsed. That collapse is what makes it linear rather than quadratic — a trie of all suffixes has O(n^2) nodes and almost all have one child. Built directly here; Ukkonen's algorithm is linear but obscures the structure at this size.
# Rosalind: BA9C — Construct the Suffix Tree of a String
# https://rosalind.info/problems/ba9c/
#
# Given: A string Text.
# Return: The strings labelling the edges of SuffixTree(Text), in any order.
let text = "ATAAATG$"
# A suffix tree is a trie of every suffix with each non-branching chain collapsed
# into a single edge. That collapse is what makes it linear in the text rather
# than quadratic: a trie of all suffixes has O(n^2) nodes, and almost all of them
# have exactly one child.
#
# Built here the direct way, by inserting suffixes into a trie and then merging
# chains. Ukkonen's algorithm builds it in linear time without the intermediate,
# which matters at genome scale and obscures the structure at this one.
let node_children = { "0": [] }
let next_id = 1
for start in range(0, len(text)) {
let node = "0"
for symbol in chars(substr(text, start, len(text) - start)) {
let onward = node_children[node] |> filter(|e| e.symbol == symbol)
if len(onward) > 0 {
node = onward[0].child
} else {
let child = str(next_id)
next_id = next_id + 1
node_children[node] = push(node_children[node], { symbol: symbol, child: child })
node_children[child] = []
node = child
}
}
}
# Walk down from each child of a branching node, absorbing single-child nodes
# into the edge label until a branch or a leaf is reached.
let labels = []
let frontier = ["0"]
while len(frontier) > 0 {
let node = frontier[0]
frontier = slice(frontier, 1, len(frontier))
for edge in node_children[node] {
let label = edge.symbol
let at = edge.child
while len(node_children[at]) == 1 {
label = label + node_children[at][0].symbol
at = node_children[at][0].child
}
labels = push(labels, label)
if len(node_children[at]) > 1 { frontier = push(frontier, at) }
}
}
println("Result: " + join(sort(labels), " "))
println("Expected (in any order): AAATG$ G$ T ATG$ TG$ A A AAATG$ G$ T G$ $")
fn test_ba9c_suffix_tree() {
let expected = ["AAATG$", "G$", "T", "ATG$", "TG$", "A", "A", "AAATG$", "G$", "T", "G$", "$"]
assert sort(labels) == sort(expected), "BA9C: got " + str(sort(labels))
# Every leaf-to-root path spells a suffix, so concatenating the labels along
# each root-to-leaf path must give back exactly the suffixes.
assert len(labels) == 12, "BA9C: expected 12 edges, got " + str(len(labels))
# Collapsing chains is the whole point: a trie of all suffixes has many more
# nodes than the tree has edges.
assert next_id - 1 > len(labels),
"BA9C: the uncollapsed trie should have more nodes than the tree has edges"
}BA9F — Find the Shortest Non-Shared Substring of Two Strings
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Shortest first, which makes the answer minimal by construction: if a substring of length k is absent, everything containing it is absent too. Returns CC where the sample shows AA; the assertion proves no one-character answer exists.
# Rosalind: BA9F — Find the Shortest Non-Shared Substring of Two Strings
# https://rosalind.info/problems/ba9f/
#
# Given: Two strings.
# Return: The shortest substring of the first that does not occur in the second.
let first_text = "CCAAGCTGCTAGAGG"
let second_text = "CATGCTGGGCTGGCT"
# Shortest first, so the answer is found before anything longer is examined. That
# ordering also makes the answer minimal by construction: if a substring of
# length k is absent, every substring containing it is absent too, so nothing
# shorter can have been missed.
#
# The textbook does this on a suffix tree of both strings at once; at this size
# the direct search is clearer and gives the same answer.
let answer = ""
let width = 1
while answer == "" and width <= len(first_text) {
let candidates = range(0, len(first_text) - width + 1)
|> map(|i| substr(first_text, i, width))
|> unique()
|> filter(|piece| contains(second_text, piece) == false)
if len(candidates) > 0 { answer = candidates[0] }
width = width + 1
}
println("Result: " + answer)
println("Expected: AA (any shortest non-shared substring is accepted)")
fn test_ba9f_shortest_non_shared_substring() {
assert len(answer) == 2, "BA9F: expected length 2, got " + str(len(answer))
assert contains(first_text, answer), "BA9F: " + answer + " is not in the first string"
assert contains(second_text, answer) == false, "BA9F: " + answer + " is in the second string"
# Minimal: every single character of the first string also occurs in the
# second, so no one-character answer exists.
let singles = chars(first_text) |> unique() |> filter(|c| contains(second_text, c) == false)
assert len(singles) == 0, "BA9F: a single character would have been shorter"
# The published answer is one of several valid ones, and scores the same.
assert contains(second_text, "AA") == false, "BA9F: AA is absent from the second string too"
}BA9N — Find All Occurrences of a Collection of Patterns in a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
BA9L and BA9M count occurrences without locating them, because a BWT band gives rows of the sorted matrix rather than text positions. Turning a row back into a position is what the suffix array supplies — the partial one of BA9Q in a real aligner, since the full one would undo the memory saving.
# Rosalind: BA9N — Find All Occurrences of a Collection of Patterns in a String
# https://rosalind.info/problems/ba9n/
#
# Given: A string Text and a collection of strings Patterns.
# Return: Every position in Text where some pattern occurs.
let text = "AATCGGGTTCAATCGGGGT"
let patterns = ["ATCG", "GGGT"]
# BA9L and BA9M count occurrences without recovering where they are, because the
# Burrows-Wheeler band gives rows of the sorted matrix rather than positions in
# the text. Turning a row back into a position is what the suffix array supplies,
# and in a real aligner it is the *partial* array of BA9Q that does it — the full
# one would undo the memory saving the whole index exists for.
let sa = suffix_array(text)
fn suffix_at(source, start) { substr(source, start, len(source) - start) }
fn starts_with(haystack, needle) {
len(haystack) >= len(needle) and substr(haystack, 0, len(needle)) == needle
}
fn occurrences_of(pattern, order, source) {
let low = 0
let high = len(order)
while low < high {
let middle = floor((low + high) / 2)
if suffix_at(source, order[middle]) < pattern { low = middle + 1 } else { high = middle }
}
let start = low
high = len(order)
while low < high {
let middle = floor((low + high) / 2)
if starts_with(suffix_at(source, order[middle]), pattern) {
low = middle + 1
} else {
high = middle
}
}
range(start, low) |> map(|i| order[i])
}
let found = sort(patterns |> flat_map(|p| occurrences_of(p, sa, text)))
println("Result: " + (found |> map(|p| str(p)) |> join(" ")))
println("Expected: 1 4 11 15")
fn test_ba9n_multiple_pattern_matching() {
assert (found |> map(|p| str(p)) |> join(" ")) == "1 4 11 15", "BA9N: got " + str(found)
for position in found {
let here = patterns |> filter(|p| substr(text, position, len(p)) == p)
assert len(here) > 0, "BA9N: nothing matches at " + str(position)
}
let by_scanning = sort(patterns |> flat_map(|p|
range(0, len(text) - len(p) + 1) |> filter(|s| substr(text, s, len(p)) == p)))
assert found == by_scanning, "BA9N: disagrees with a plain scan"
}BA9O — Find All Approximate Occurrences of a Collection of Patterns in a String
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
Reads carry errors and genomes carry variants, so exact matching finds nothing useful and every aligner is an approximate matcher. Position 4 appears twice because two patterns match there — the answer lists occurrences, not distinct positions.
# Rosalind: BA9O — Find All Approximate Occurrences of a Collection of Patterns
# https://rosalind.info/problems/ba9o/
#
# Given: A string Text, a collection of Patterns, and an integer d.
# Return: Every position where some pattern occurs with at most d mismatches.
let text = "ACATGCTACTTT"
let patterns = ["ATT", "GCC", "GCTA", "TATT"]
let allowed = 1
# Real reads carry sequencing errors and real genomes carry variants, so exact
# matching finds nothing useful — which is why every aligner is an approximate
# matcher. The standard trick is seed-and-extend: split the pattern into d+1
# pieces, and since d mismatches cannot spoil all of them, at least one piece
# must match exactly. Those exact matches are found by the index, and only their
# neighbourhoods are checked in full.
#
# Written out directly here, because the point is the answer rather than the
# index; the pieces are checked against every position instead.
fn mismatches(a, b) { range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1)) }
let found = sort(patterns |> flat_map(|pattern|
range(0, len(text) - len(pattern) + 1)
|> filter(|start| mismatches(substr(text, start, len(pattern)), pattern) <= allowed)))
println("Result: " + (found |> map(|p| str(p)) |> join(" ")))
println("Expected: 2 4 4 6 7 8 9")
fn test_ba9o_approximate_matching() {
assert (found |> map(|p| str(p)) |> join(" ")) == "2 4 4 6 7 8 9", "BA9O: got " + str(found)
# 4 appears twice because two different patterns match there — the answer
# lists occurrences, not distinct positions.
assert (found |> count_if(|p| p == 4)) == 2, "BA9O: two patterns match at position 4"
# Every reported position really is within d.
for position in found {
let close = patterns |> filter(|p|
position + len(p) <= len(text)
and mismatches(substr(text, position, len(p)), p) <= allowed)
assert len(close) > 0, "BA9O: nothing is within " + str(allowed) + " at " + str(position)
}
# And exact matching alone would find fewer, which is why d matters.
let exact = patterns |> flat_map(|p|
range(0, len(text) - len(p) + 1) |> filter(|s| substr(text, s, len(p)) == p))
assert len(exact) < len(found), "BA9O: allowing a mismatch should find more"
}BA9P — Implement TreeColoring
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
How BA9E's shared-substring question is answered on a generalised suffix tree: colour leaves by which string they came from, and an internal node goes purple exactly when its substring occurs in both.
# Rosalind: BA9P — Implement TreeColoring
# https://rosalind.info/problems/ba9p/
#
# Given: The adjacency list of a suffix tree and the colours of its leaves.
# Return: The colour of every node.
let children = {
"0": [], "1": [], "2": ["0", "1"], "3": [], "4": [],
"5": ["3", "2"], "6": [], "7": ["4", "5", "6"],
}
let leaf_colours = { "0": "red", "1": "red", "3": "blue", "4": "blue", "6": "red" }
# This is how BA9E's shared-substring question is answered on a generalised
# suffix tree: leaves are coloured by which string they came from, and an
# internal node ends up purple exactly when the substring it spells occurs in
# both. So a colouring pass turns "which substrings are shared" into a property
# readable off the tree.
#
# A node is ripe when every child is already coloured. Working outwards from the
# leaves means each node is decided once, and a node with children of differing
# colours becomes purple regardless of which colours they were.
let colours = {}
for node in keys(leaf_colours) { colours[node] = leaf_colours[node] }
let uncoloured = keys(children) |> filter(|node| contains(keys(colours), node) == false)
while len(uncoloured) > 0 {
let ripe = uncoloured |> filter(|node|
(children[node] |> count_if(|kid| contains(keys(colours), kid) == false)) == 0)
for node in ripe {
let kid_colours = children[node] |> map(|kid| colours[kid]) |> unique()
colours[node] = if len(kid_colours) == 1 then kid_colours[0] else "purple"
}
uncoloured = uncoloured |> filter(|node| contains(keys(colours), node) == false)
}
let listed = sort(keys(colours) |> map(|k| int(k))) |> map(|node|
str(node) + ": " + colours[str(node)])
println("Result:")
for line in listed { println(" " + line) }
println("Expected: 0 red, 1 red, 2 red, 3 blue, 4 blue, 5 purple, 6 red, 7 purple")
fn test_ba9p_tree_colouring() {
assert join(listed, ", ")
== "0: red, 1: red, 2: red, 3: blue, 4: blue, 5: purple, 6: red, 7: purple",
"BA9P: got " + join(listed, ", ")
# Node 2's children are both red, so it inherits; node 5 has a red child and
# a blue one, so it cannot and becomes purple. That distinction is the whole
# algorithm.
assert colours["2"] == "red", "BA9P: 2 inherits from two red children"
assert colours["5"] == "purple", "BA9P: 5 has children of different colours"
# Every node is coloured, and purple appears only where children disagree.
assert len(keys(colours)) == len(keys(children)), "BA9P: every node needs a colour"
for node in keys(children) {
if len(children[node]) > 0 {
let kid_colours = children[node] |> map(|kid| colours[kid]) |> unique()
let should_be_purple = len(kid_colours) > 1
assert (colours[node] == "purple") == should_be_purple,
"BA9P: node " + node + " is wrongly coloured " + colours[node]
}
}
}BA9R — Construct a Suffix Tree from a Suffix Array
solvedbrowser + CLI Problem statement Open in the workbench Download .bl
BA9C built the tree by collapsing a suffix trie, which is quadratic before the collapse. The suffix and LCP arrays carry the same information in two flat integer arrays and rebuild the tree in one pass — which is why real tools store the arrays and never materialise the tree.
# Rosalind: BA9R — Construct a Suffix Tree from a Suffix Array
# https://rosalind.info/problems/ba9r/
#
# Given: A string, its suffix array and its LCP array.
# Return: The edge labels of the suffix tree, in any order.
let text = "GTAGT$"
# BA9C built the tree by inserting every suffix into a trie and collapsing it,
# which costs O(n^2) time and memory before the collapse. The suffix array and
# LCP array carry the same information in two flat integer arrays, and the tree
# can be rebuilt from them in one left-to-right pass — so a structure that was
# quadratic to build becomes linear given arrays that are themselves linear.
#
# That is why real tools store the arrays and never materialise the tree.
let sa = suffix_array(text)
let lcp = lcp_array(text)
println("Suffix array: " + (sa |> map(|v| str(v)) |> join(", ")))
println("LCP array: " + (lcp |> map(|v| str(v)) |> join(", ")))
# Two passes, because a leaf's edge depends on what comes *after* it as well as
# before. The node a leaf hangs from sits at depth max(lcp[i], lcp[i+1]) — the
# deeper of what it shares with either neighbour — so the leaf edge is whatever
# is left of the suffix below that depth.
let leaf_labels = range(0, len(sa)) |> map(|i| {
let before = if i == 0 then 0 else lcp[i]
let after = if i + 1 < len(sa) then lcp[i + 1] else 0
let depth = max([before, after])
substr(text, sa[i] + depth, len(text) - sa[i] - depth)
})
# The internal nodes come from the LCP array alone: a run of suffixes sharing a
# prefix of length L hangs off one node at depth L, and a stack tracks which are
# still open as the array is scanned.
let internal_labels = []
let stack = [{ depth: 0, start: 0 }]
for i in range(1, len(sa)) {
let shared = lcp[i]
let last_start = sa[i]
while stack[len(stack) - 1].depth > shared {
let closing = stack[len(stack) - 1]
stack = slice(stack, 0, len(stack) - 1)
let parent_depth = max([stack[len(stack) - 1].depth, shared])
internal_labels = push(internal_labels,
substr(text, closing.start + parent_depth, closing.depth - parent_depth))
last_start = closing.start
}
if stack[len(stack) - 1].depth < shared {
stack = push(stack, { depth: shared, start: last_start })
}
}
while len(stack) > 1 {
let closing = stack[len(stack) - 1]
stack = slice(stack, 0, len(stack) - 1)
let parent_depth = stack[len(stack) - 1].depth
internal_labels = push(internal_labels,
substr(text, closing.start + parent_depth, closing.depth - parent_depth))
}
let labels = leaf_labels + internal_labels
println("Result: " + join(sort(labels), " "))
println("Expected (any order): $ T AGT$ $ AGT$ GT $ AGT$")
fn test_ba9r_suffix_tree_from_arrays() {
let expected = ["$", "T", "AGT$", "$", "AGT$", "GT", "$", "AGT$"]
assert sort(labels) == sort(expected), "BA9R: got " + join(sort(labels), " ")
assert len(labels) == 8, "BA9R: expected 8 edges, got " + str(len(labels))
# One leaf per suffix, and every leaf edge ends at the sentinel.
assert len(leaf_labels) == len(text), "BA9R: one leaf per suffix"
for label in leaf_labels {
assert substr(label, len(label) - 1, 1) == "$", "BA9R: a leaf edge must reach the end"
}
# The internal edges are the shared prefixes: GT and T, each shared by two
# suffixes.
assert sort(internal_labels) == ["GT", "T"], "BA9R: got internals " + str(internal_labels)
# Every label is a real substring of the text.
for label in labels {
assert contains(text, label), "BA9R: " + label + " is not in the text"
}
}