Every problem

DNA — Counting DNA Nucleotides

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A warm-up for the site's input format rather than a biological question — but the four counts it produces are exactly what GC content is computed from, and a lopsided base composition is the first thing a QC tool flags as contamination or adapter read-through.

# Rosalind: DNA — Counting DNA Nucleotides
# https://rosalind.info/problems/dna/
#
# Given: A DNA string s of length at most 1000 nt.
# Return: Four integers counting A, C, G, T.

let s = dna"AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC"

let counts = base_counts(s)
let result = str(counts.A) + " " + str(counts.C) + " " + str(counts.G) + " " + str(counts.T)

println("Result:   " + result)
println("Expected: 20 12 17 21")

fn test_dna_nucleotide_counts() {
    assert result == "20 12 17 21", "DNA: got '" + result + "'"
}

RNA — Transcribing DNA into RNA

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

One character substituted for another, which is a fair summary of transcription only if you ignore splicing, capping and polyadenylation. SPLC covers the part this leaves out.

# Rosalind: RNA — Transcribing DNA into RNA
# https://rosalind.info/problems/rna/
#
# Given: A DNA string t.
# Return: The RNA string u, with every T replaced by U.

let t = dna"GATGGAACTTGACTACGTAAATT"
let u = transcribe(t)

println("Result:   " ++ str(u))
println("Expected: GAUGGAACUUGACUACGUAAAUU")

fn test_rna_transcription() {
    assert str(u) == "GAUGGAACUUGACUACGUAAAUU", "RNA: got " ++ str(u)
}

REVC — Complementing a Strand of DNA

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The most-used operation in the whole track. A read carries no record of which strand it came from, so a motif present on one strand appears as its reverse complement on the other — which is why REVP, BA6E and GASM all have to search both.

# Rosalind: REVC — Complementing a Strand of DNA
# https://rosalind.info/problems/revc/
#
# Given: A DNA string s of length at most 1000 bp.
# Return: The reverse complement of s.

let s = dna"AAAACCCGGT"
let rc = reverse_complement(s)

println("Result:   " ++ str(rc))
println("Expected: ACCGGGTTTT")

fn test_revc_reverse_complement() {
    assert str(rc) == "ACCGGGTTTT", "REVC: got " ++ str(rc)
}

HAMM — Counting Point Mutations

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Counts differing positions, which undercounts mutations: a site can change and change back, leaving no trace. RSUB shows exactly that happening, and it is why distance models correct for unseen substitutions rather than using raw counts.

# Rosalind: HAMM — Counting Point Mutations
# https://rosalind.info/problems/hamm/
#
# Given: Two DNA strings s and t of equal length.
# Return: The Hamming distance between them.

let s = dna"GAGCCTACTAACGGGAT"
let t = dna"CATCGTAATGACGGCCT"

let distance = hamming_distance(s, t)

println("Result:   " + str(distance))
println("Expected: 7")

fn test_hamm_point_mutations() {
    assert distance == 7, "HAMM: got " + str(distance)
}

FIB — Rabbits and Recurrence Relations

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A recurrence-relation exercise wearing rabbits, not population biology — nothing here dies, competes or mutates. WFMD and EBIN are where real population dynamics appear.

# Rosalind: FIB — Rabbits and Recurrence Relations
# https://rosalind.info/problems/fib/
#
# Given: Positive integers n <= 40 and k <= 5.
# Return: The total rabbit pairs after n months, when every pair of
# reproduction-age rabbits produces a litter of k pairs each month.

let n = 5
let k = 3

fn rabbit_pairs(months, litter) {
    let previous = 1
    let current = 1
    let month = 3
    while month <= months {
        let next = current + previous * litter
        previous = current
        current = next
        month = month + 1
    }
    if months <= 2 then 1 else current
}

let result = rabbit_pairs(n, k)

println("Result:   " + str(result))
println("Expected: 19")

fn test_fib_rabbit_pairs() {
    assert result == 19, "FIB: got " + str(result)
}

GC — Computing GC Content

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

GC content ranges from about 16 to 75 percent across bacteria, and a difference under 5 percent is taken as evidence of the same species — which is why a fragment's GC content alone narrows down where it came from. It also biases sequencing: GC-rich regions amplify poorly and end up under-covered.

# Rosalind: GC — Computing GC Content
# https://rosalind.info/problems/gc/
#
# Given: At most 10 DNA strings in FASTA format.
# Return: The ID of the string with the highest GC content, and that content.

let records = [
    { id: "Rosalind_6404", seq: dna"CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAATAATTCTGAGG" },
    { id: "Rosalind_5959", seq: dna"CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTTGTCAGCAGACACGC" },
    { id: "Rosalind_0808", seq: dna"CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT" }
]

let scored = records |> map(|r| { id: r.id, gc: gc_content(r.seq) * 100.0 })
let best = scored |> sort_by(|r| r.gc) |> reverse() |> first()

println("Result:   " + best.id + " " + str(best.gc))
println("Expected: Rosalind_0808 60.919540")

fn test_gc_highest_content() {
    assert best.id == "Rosalind_0808", "GC: picked " + best.id
    assert round(best.gc, 6) == 60.91954, "GC: got " + str(best.gc)
}

PROT — Translating RNA into Protein

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Sixty-four codons encode twenty amino acids and a stop, so the code is degenerate and translation discards information. MRNA counts exactly how much.

# Rosalind: PROT — Translating RNA into Protein
# https://rosalind.info/problems/prot/
#
# Given: An RNA string s corresponding to a strand of mRNA.
# Return: The protein string encoded by s.

let s = rna"AUGGCCAUGGCGCCCAGAACUGAGAUCAAUAGUACCCGUAUUAACGGGUGA"

# translate() stops at the first stop codon, which is what the problem asks for.
let peptide = translate(s)

println("Result:   " ++ str(peptide))
println("Expected: MAMAPRTEINSTRING")

fn test_prot_translation() {
    assert str(peptide) == "MAMAPRTEINSTRING", "PROT: got " ++ str(peptide)
}

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 + "'"
}

FIBD — Mortal Fibonacci Rabbits

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

FIB with a death clock added. Still a recurrence exercise rather than population biology, and useful mainly for showing that the state has to widen from one number to an age profile.

# Rosalind: FIBD — Mortal Fibonacci Rabbits
# https://rosalind.info/problems/fibd/
#
# Given: Positive integers n <= 100 and m <= 20.
# Return: The number of rabbit pairs alive after n months, when every pair
# lives for exactly m months.

let n = 6
let m = 3

# ages[i] = pairs currently i months old. Each month the adults breed and the
# oldest cohort dies.
fn mortal_pairs(months, lifespan) {
    let ages = range(0, lifespan) |> map(|i| if i == 0 then 1 else 0)
    let month = 1
    while month < months {
        let newborns = sum(drop(ages, 1))
        ages = concat([newborns], take(ages, lifespan - 1))
        month = month + 1
    }
    sum(ages)
}

let result = mortal_pairs(n, m)

println("Result:   " + str(result))
println("Expected: 4")

fn test_fibd_mortal_rabbits() {
    assert result == 4, "FIBD: got " + str(result)
}

IPRB — Mendel's First Law

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Mendel's first law as a probability question. The whole content is that two alleles segregate independently, which is what makes the offspring distribution a product rather than something requiring simulation.

# Rosalind: IPRB — Mendel's First Law
# https://rosalind.info/problems/iprb/
#
# Given: Positive integers k (homozygous dominant), m (heterozygous) and
# n (homozygous recessive) organisms.
# Return: The probability that two randomly chosen mating organisms produce an
# individual with a dominant allele.

let k = 2.0
let m = 2.0
let n = 2.0
let total = k + m + n

# Work out the recessive probability and subtract: only mm, mn and nn pairings
# can produce a recessive offspring.
# A continued expression keeps its operator at the end of the line.
let mm = (m / total) * ((m - 1.0) / (total - 1.0)) * 0.25
let mn = (m / total) * (n / (total - 1.0)) * 0.5
let nm = (n / total) * (m / (total - 1.0)) * 0.5
let nn = (n / total) * ((n - 1.0) / (total - 1.0))
let p_recessive = mm + mn + nm + nn

let result = 1.0 - p_recessive

println("Result:   " + str(round(result, 5)))
println("Expected: 0.78333")

fn test_iprb_dominant_probability() {
    assert round(result, 5) == 0.78333, "IPRB: got " + str(result)
}

IEV — Calculating Expected Offspring

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Expectation is linear whether or not the events are independent, so the answer is a weighted sum with no interaction between genotypes to model. That property is why expected values are reached for far more often than full distributions.

# Rosalind: IEV — Calculating Expected Offspring
# https://rosalind.info/problems/iev/
#
# Given: Six integers, the number of couples with each genotype pairing.
# Return: The expected number of offspring displaying the dominant phenotype,
# assuming every couple has exactly two offspring.

let couples = [1, 0, 0, 1, 0, 1]

# Probability a single offspring shows the dominant phenotype, per pairing:
# AA-AA, AA-Aa, AA-aa are certain; Aa-Aa is 3/4; Aa-aa is 1/2; aa-aa is 0.
let dominant_probability = [1.0, 1.0, 1.0, 0.75, 0.5, 0.0]
let offspring_per_couple = 2.0

let result = range(0, 6)
    |> map(|i| float(couples[i]) * dominant_probability[i] * offspring_per_couple)
    |> sum()

println("Result:   " + str(result))
println("Expected: 3.5")

fn test_iev_expected_offspring() {
    assert result == 3.5, "IEV: got " + str(result)
}

MRNA — Inferring mRNA from Protein

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The reverse of PROT, and it cannot be done — the code is degenerate, so a protein maps back to many mRNAs. Leucine, serine and arginine take six codons each. Counting them needs a modulus because the total outgrows any integer.

# Rosalind: MRNA — Inferring mRNA from Protein
# https://rosalind.info/problems/mrna/
#
# Given: A protein string of length at most 1000 aa.
# Return: The number of RNA strings that could have produced it, modulo
# 1,000,000 — remembering the stop codon.

let protein_string = "MA"

# Codons per amino acid in the standard genetic code.
let codon_counts = {
    A: 4, C: 2, D: 2, E: 2, F: 2, G: 4, H: 2, I: 3, K: 2, L: 6,
    M: 1, N: 2, P: 4, Q: 2, R: 6, S: 6, T: 4, V: 4, W: 1, Y: 2
}
let stop_codons = 3
let modulus = 1000000

let result = range(0, len(protein_string))
    |> map(|i| codon_counts[substr(protein_string, i, 1)])
    |> reduce(|acc, n| (acc * n) % modulus, stop_codons)

println("Result:   " + str(result))
println("Expected: 12")

fn test_mrna_possible_rna_strings() {
    assert result == 12, "MRNA: got " + str(result)
}

PRTM — Calculating Protein Mass

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Monoisotopic masses, which is what a high-resolution instrument reports when it can resolve individual isotope peaks; low-resolution and high-mass measurements give average masses instead, differing by roughly 0.06 percent. Everything in the mass-spectrometry chapter rests on this sum.

# Rosalind: PRTM — Calculating Protein Mass
# https://rosalind.info/problems/prtm/
#
# Given: A protein string P of length at most 1000 aa.
# Return: The total weight of P, using the monoisotopic mass table.

let protein_string = "SKADYEK"

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}

let result = range(0, len(protein_string))
    |> map(|i| monoisotopic[substr(protein_string, i, 1)])
    |> sum()

println("Result:   " + str(round(result, 3)))
println("Expected: 821.392")

fn test_prtm_protein_mass() {
    assert round(result, 3) == 821.392, "PRTM: got " + str(result)
}

PERM — Enumerating Gene Orders

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Permutations as gene orders. n! grows fast enough that enumeration stops being possible almost immediately, which is the premise SIGN, REAR and SORT build on — there, gene order is the data and the question is how few reversals separate two of them.

# Rosalind: PERM — Enumerating Gene Orders
# https://rosalind.info/problems/perm/
#
# Given: A positive integer n <= 7.
# Return: The total number of permutations of length n, followed by a list of
# all such permutations.

let n = 3

fn permutations_of(items) {
    if len(items) <= 1 {
        [items]
    } else {
        range(0, len(items)) |> flat_map(|i| {
            let chosen = items[i]
            let rest = concat(take(items, i), drop(items, i + 1))
            permutations_of(rest) |> map(|p| concat([chosen], p))
        })
    }
}

let orders = permutations_of(range(1, n + 1))

println("Result:   " + str(len(orders)))
println("Expected: 6")
orders |> each(|p| println("  " + (p |> map(|x| str(x)) |> join(" "))))

fn test_perm_permutation_count() {
    assert len(orders) == 6, "PERM: got " + str(len(orders))
    assert len(unique(orders |> map(|p| join(map(p, |x| str(x)), " ")))) == 6, "PERM: duplicates"
}

PPER — Partial Permutations

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Ordered selections rather than full orderings, taken modulo a million because the count outgrows the answer format long before it outgrows the biology. Combinatorial groundwork for the counting problems, not a biological question in itself.

# Rosalind: PPER — Partial Permutations
# https://rosalind.info/problems/pper/
#
# Given: Positive integers n and k with 100 >= n > 0 and 10 >= k > 0.
# Return: The number of partial permutations P(n, k), modulo 1,000,000.

let n = 21
let k = 7
let modulus = 1000000

# P(n, k) = n * (n-1) * ... * (n-k+1), reduced at every step to stay small.
let result = range(0, k) |> reduce(|acc, i| (acc * (n - i)) % modulus, 1)

println("Result:   " + str(result))
println("Expected: 51200")

fn test_pper_partial_permutations() {
    assert result == 51200, "PPER: got " + str(result)
}

EDIT — Edit Distance

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The crude version of sequence comparison, and the pack refines it twice: GLOB stops treating every substitution as equally costly, and GAFF stops treating a gap of k bases as k separate events. Useful mainly as the baseline those two improve on.

# Rosalind: EDIT — Edit Distance
# https://rosalind.info/problems/edit/
#
# Given: Two protein strings s and t.
# Return: The edit distance between them.

let s = "PLEASANTLY"
let t = "MEANLY"

let distance = edit_distance(s, t)

println("Result:   " + str(distance))
println("Expected: 5")

fn test_edit_distance() {
    assert distance == 5, "EDIT: got " + str(distance)
}

GLOB — Global Alignment with Scoring Matrix

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

BLOSUM62 is built from substitution frequencies actually observed in aligned protein blocks, so swapping a residue for a chemically similar one costs little and swapping it for an unrelated one costs a lot. That is the information EDIT throws away by charging one for every mismatch.

# Rosalind: GLOB — Global Alignment with Scoring Matrix
# https://rosalind.info/problems/glob/
#
# Given: Two protein strings.
# Return: The maximum global alignment score using BLOSUM62 and a constant
# gap penalty of 5.

let s = "PLEASANTLY"
let t = "MEANLY"
let gap = -5.0

let blosum = score_matrix("blosum62")
let residues = blosum.row_names
let width = blosum.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# Needleman-Wunsch with a linear gap penalty. Rows walk s, columns walk t.
fn global_score(a, b, mat, names, cols, gap_penalty) {
    let previous = range(0, len(b) + 1) |> map(|j| float(j) * gap_penalty)
    let i = 1
    while i <= len(a) {
        let current = [float(i) * gap_penalty]
        let ai = substr(a, i - 1, 1)
        let j = 1
        while j <= len(b) {
            let diagonal = previous[j - 1] + substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let up = previous[j] + gap_penalty
            let left = current[j - 1] + gap_penalty
            current = push(current, max([diagonal, up, left]))
            j = j + 1
        }
        previous = current
        i = i + 1
    }
    previous[len(b)]
}

let result = global_score(s, t, blosum, residues, width, gap)

println("Result:   " + str(int(result)))
println("Expected: 8")

fn test_glob_blosum62_global_alignment() {
    assert int(result) == 8, "GLOB: got " + str(result)
}

LOCA — Local Alignment with Scoring Matrix

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Two proteins may share a single conserved domain inside otherwise unrelated sequence. A global alignment has to align the unrelated parts too and buries the signal; local alignment reports only the stretch that matches, which is what makes it the basis of database search.

# Rosalind: LOCA — Local Alignment with Scoring Matrix
# https://rosalind.info/problems/loca/
#
# Given: Two protein strings.
# Return: The maximum local alignment score using PAM250 and a constant gap
# penalty of 5.

let s = "MEANLYPRTEINSTRING"
let t = "PLEASANTLYEINSTEIN"
let gap = -5.0

let pam = score_matrix("pam250")
let residues = pam.row_names
let width = pam.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# Smith-Waterman: identical recurrence to the global case except that a cell
# never drops below zero, and the answer is the best cell anywhere.
fn local_score(a, b, mat, names, cols, gap_penalty) {
    let previous = range(0, len(b) + 1) |> map(|j| 0.0)
    let best = 0.0
    let i = 1
    while i <= len(a) {
        let current = [0.0]
        let ai = substr(a, i - 1, 1)
        let j = 1
        while j <= len(b) {
            let diagonal = previous[j - 1] + substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let cell = max([0.0, diagonal, previous[j] + gap_penalty, current[j - 1] + gap_penalty])
            current = push(current, cell)
            if cell > best then best = cell
            j = j + 1
        }
        previous = current
        i = i + 1
    }
    best
}

let result = local_score(s, t, pam, residues, width, gap)

println("Result:   " + str(int(result)))
println("Expected: 23")

fn test_loca_pam250_local_alignment() {
    assert int(result) == 23, "LOCA: got " + str(result)
}

GAFF — Global Alignment with Scoring Matrix and Affine Gap Penalty

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

One insertion of ten bases is a single mutational event, not ten. Charging per base makes an aligner scatter several short gaps where one long gap is right, so opening a gap costs more than extending it.

# Rosalind: GAFF — Global Alignment with Scoring Matrix and Affine Gap Penalty
# https://rosalind.info/problems/gaff/
#
# Given: Two protein strings.
# Return: The maximum global alignment score using BLOSUM62, a gap opening
# penalty of 11 and a gap extension penalty of 1.

let s = "PRTEINS"
let t = "PRTWPSEIN"
let gap_open = -11.0
let gap_extend = -1.0

let blosum = score_matrix("blosum62")
let residues = blosum.row_names
let width = blosum.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# Gotoh's three matrices, held one row at a time:
#   M — s[i] aligned to t[j]
#   X — a gap in t (consuming s)
#   Y — a gap in s (consuming t)
# Opening a gap costs gap_open, each further residue costs gap_extend.
fn affine_score(a, b, mat, names, cols, open_penalty, extend_penalty) {
    let n = len(b)
    let very_low = -1000000.0

    let prev_m = concat([0.0], range(1, n + 1) |> map(|_| very_low))
    let prev_x = concat([very_low], range(1, n + 1) |> map(|_| very_low))
    let prev_y = concat([very_low], range(1, n + 1) |> map(|j| open_penalty + float(j - 1) * extend_penalty))

    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        let row_m = [very_low]
        let row_x = [open_penalty + float(i - 1) * extend_penalty]
        let row_y = [very_low]

        let j = 1
        while j <= n {
            let sub = substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let best_prev = max([prev_m[j - 1], prev_x[j - 1], prev_y[j - 1]])
            row_m = push(row_m, best_prev + sub)
            row_x = push(row_x, max([prev_m[j] + open_penalty, prev_x[j] + extend_penalty]))
            row_y = push(row_y, max([row_m[j - 1] + open_penalty, row_y[j - 1] + extend_penalty]))
            j = j + 1
        }

        prev_m = row_m
        prev_x = row_x
        prev_y = row_y
        i = i + 1
    }

    max([prev_m[n], prev_x[n], prev_y[n]])
}

let result = affine_score(s, t, blosum, residues, width, gap_open, gap_extend)

println("Result:   " + str(int(result)))
println("Expected: 8")

fn test_gaff_affine_global_alignment() {
    assert int(result) == 8, "GAFF: got " + str(result)
}

OAP — Overlap Alignment

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Brute-force DP over lists; ~9 s. A flat-buffer implementation would be the natural follow-up.

# Rosalind: OAP — Overlap Alignment
# https://rosalind.info/problems/oap/
#
# Given: Two DNA strings s and t.
# Return: The score of an optimal overlap alignment of s and t, where a suffix
# of s is aligned against a prefix of t. Match +1, mismatch and gap -2.

let s = "CTAAGGGATTCCGGTAATTAGACAG"
let t = "ATAGACCATATGTCAGTGACTGTGTAA"
let match_score = 1.0
let penalty = -2.0

# An overlap alignment is a global alignment with two changes: starting
# anywhere in s is free (first column zero), and ending anywhere in t is free
# (the answer is the best value in the last row).
fn overlap_score(a, b, match_value, mismatch_value) {
    let n = len(b)
    let previous = range(0, n + 1) |> map(|j| float(j) * mismatch_value)

    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        # Free start in s: no penalty accumulated down the first column.
        let current = [0.0]
        let j = 1
        while j <= n {
            let same = ai == substr(b, j - 1, 1)
            let diagonal = previous[j - 1] + (if same then match_value else mismatch_value)
            current = push(current, max([diagonal, previous[j] + mismatch_value, current[j - 1] + mismatch_value]))
            j = j + 1
        }
        previous = current
        i = i + 1
    }

    # Free end in t: the best score anywhere along the final row.
    max(previous)
}

let result = overlap_score(s, t, match_score, penalty)

println("Result:   " + str(int(result)))
println("Expected: 1")

fn test_oap_overlap_alignment() {
    assert int(result) == 1, "OAP: got " + str(result)
}

SSET — Counting Subsets

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Two to the n, which is the point. Any method that examines every subset of a set of sites stops being usable in the low twenties, and that is the wall the search problems elsewhere in this pack are built to avoid.

# Rosalind: SSET — Counting Subsets
# https://rosalind.info/problems/sset/
#
# Given: A positive integer n <= 1000.
# Return: The total number of subsets of {1, 2, ..., n}, modulo 1,000,000.

let n = 3
let modulus = 1000000

# 2^n, reduced at every step so the value never grows beyond the modulus.
let result = range(0, n) |> reduce(|acc, _| (acc * 2) % modulus, 1)

println("Result:   " + str(result))
println("Expected: 8")

fn test_sset_subset_count() {
    assert result == 8, "SSET: got " + str(result)
}

PMCH — Perfect Matchings and RNA Secondary Structures

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

RNA folds back on itself and pairs. This counts the ways every base could be paired, which requires the A and U counts to match and the C and G counts to match — a condition real sequences rarely satisfy, which is what MMCH addresses.

# Rosalind: PMCH — Perfect Matchings and RNA Secondary Structures
# https://rosalind.info/problems/pmch/
#
# Given: An RNA string s with the same number of A as U and of G as C.
# Return: The total number of perfect matchings of basepair edges.

let s = "AGCUAGUCAU"

fn count_base(seq, base) {
    range(0, len(seq)) |> count_if(|i| substr(seq, i, 1) == base)
}

fn factorial(n) {
    range(1, n + 1) |> reduce(|acc, i| acc * i, 1)
}

# Every A pairs with some U and every G with some C, independently, so the
# count is |A|! * |G|!.
let adenine = count_base(s, "A")
let guanine = count_base(s, "G")
let result = factorial(adenine) * factorial(guanine)

println("Result:   " + str(result))
println("Expected: 12")

fn test_pmch_perfect_matchings() {
    assert result == 12, "PMCH: got " + str(result)
}

MMCH — Maximum Matchings and RNA Secondary Structures

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

PMCH without the assumption that everything pairs. Real RNA leaves bases unpaired because the counts do not match, so the question becomes how many bonds are possible rather than how many arrangements are complete.

# Rosalind: MMCH — Maximum Matchings and RNA Secondary Structures
# https://rosalind.info/problems/mmch/
#
# Given: An RNA string s.
# Return: The total number of maximum matchings of basepair edges.

let s = "AUGCUUC"

fn count_base(seq, base) {
    range(0, len(seq)) |> count_if(|i| substr(seq, i, 1) == base)
}

# P(n, k) — the number of ways to pair k of the scarcer base with n of the
# commoner one.
fn partial_permutations(n, k) {
    range(0, k) |> reduce(|acc, i| acc * (n - i), 1)
}

fn pairing_ways(first, second) {
    partial_permutations(max([first, second]), min([first, second]))
}

# A continued expression keeps its operator at the end of the line, so these are
# bound separately rather than split across a leading `*`.
let au_ways = pairing_ways(count_base(s, "A"), count_base(s, "U"))
let gc_ways = pairing_ways(count_base(s, "G"), count_base(s, "C"))
let result = au_ways * gc_ways

println("Result:   " + str(result))
println("Expected: 6")

fn test_mmch_maximum_matchings() {
    assert result == 6, "MMCH: got " + str(result)
}

CAT — Catalan Numbers and RNA Secondary Structures

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Catalan numbers count matchings that do not cross, and non-crossing is a physical constraint rather than a mathematical convenience: crossing bonds are pseudoknots, which this model excludes. MOTZ relaxes the requirement that everything pairs, and RNAS adds wobble bonds and a minimum hairpin turn.

# Rosalind: CAT — Catalan Numbers and RNA Secondary Structures
# https://rosalind.info/problems/cat/
#
# Given: An RNA string s.
# Return: The total number of noncrossing perfect matchings of basepair edges,
# modulo 1,000,000.

let s = "AUAU"
let modulus = 1000000

fn complements(a, b) {
    (a == "A" and b == "U") or (a == "U" and b == "A") or (a == "G" and b == "C") or (a == "C" and b == "G")
}

# Noncrossing means the first base pairs with some k, splitting the rest into
# an inside and an outside that cannot interleave — so the counts multiply.
# Only even-length intervals can be matched at all.
fn noncrossing(seq, lo, hi) {
    let width = hi - lo
    if width <= 0 {
        1
    } else {
        let first = substr(seq, lo, 1)
        let total = range(lo + 1, hi)
            |> filter(|k| ((k - lo) % 2 == 1) and complements(first, substr(seq, k, 1)))
            |> map(|k| (noncrossing(seq, lo + 1, k) * noncrossing(seq, k + 1, hi)) % modulus)
            |> sum()
        total % modulus
    }
}

let result = noncrossing(s, 0, len(s))

println("Result:   " + str(result))
println("Expected: 2")

fn test_cat_noncrossing_matchings() {
    assert result == 2, "CAT: got " + str(result)
}

LIA — Independent Alleles

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Mendel's second law. Alleles at different loci assort independently, so the two-locus probability factors — which is what makes this a binomial calculation rather than a simulation.

# Rosalind: LIA — Independent Alleles
# https://rosalind.info/problems/lia/
#
# Given: Integers k <= 7 and N <= 2^k. Tom is AaBb; every organism mates with
# an AaBb partner and has two children.
# Return: The probability that at least N organisms in generation k are AaBb.

let k = 2
let n = 1

# Mendel's second law makes the two loci independent, so each child is AaBb
# with probability 1/4 regardless of generation. Generation k holds 2^k
# organisms, so the count is Binomial(2^k, 1/4).
let population = int(pow(2.0, float(k)))
let p = 0.25

fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) }
fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) }

# P(at least N) = 1 - P(fewer than N).
let below = range(0, n)
    |> map(|i| float(choose(population, i)) * pow(p, float(i)) * pow(1.0 - p, float(population - i)))
    |> sum()

let result = 1.0 - below

println("Result:   " + str(round(result, 3)))
println("Expected: 0.684")

fn test_lia_independent_alleles() {
    assert round(result, 3) == 0.684, "LIA: got " + str(result)
}

SEXL — Sex-Linked Inheritance

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Males carry one X, so a recessive allele on it is expressed with no second copy to mask it. The proportion of affected males is therefore the allele frequency itself, which is why X-linked recessive conditions appear far more often in males than females.

# Rosalind: SEXL — Sex-Linked Inheritance
# https://rosalind.info/problems/sexl/
#
# Given: An array A of allele frequencies for a gene on the X chromosome.
# Return: For each frequency, the probability that a randomly selected female
# is a carrier of the recessive allele.

let freqs = [0.1, 0.5, 0.8]

# Males carry one X, so the recessive allele frequency is the given value. A
# female is heterozygous with probability 2q(1-q).
let results = freqs |> map(|q| 2.0 * q * (1.0 - q))
let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0.18 0.5 0.32")

fn test_sexl_carrier_probabilities() {
    assert round(results[0], 3) == 0.18, "SEXL[0]: got " + str(results[0])
    assert round(results[1], 3) == 0.5, "SEXL[1]: got " + str(results[1])
    assert round(results[2], 3) == 0.32, "SEXL[2]: got " + str(results[2])
}

AFRQ — Counting Disease Carriers

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Hardy-Weinberg run backwards: the disease frequency gives the allele frequency, and the carrier frequency 2pq follows. For a rare recessive allele carriers vastly outnumber sufferers, which is the counterintuitive result the arithmetic exists to make obvious.

# Rosalind: AFRQ — Counting Disease Carriers
# https://rosalind.info/problems/afrq/
#
# Given: An array A giving, for each factor, the proportion of the population
# that is homozygous recessive.
# Return: The proportion carrying at least one copy of the recessive allele,
# assuming Hardy-Weinberg equilibrium.

let homozygous_recessive = [0.1, 0.25, 0.5]

# Under Hardy-Weinberg, q^2 is given, so q = sqrt(q^2) and the carriers plus
# the affected are everyone who is not homozygous dominant: 1 - (1 - q)^2.
let results = homozygous_recessive |> map(|q2| {
    let q = sqrt(q2)
    1.0 - (1.0 - q) * (1.0 - q)
})
let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0.532 0.75 0.914")

fn test_afrq_carrier_proportions() {
    assert round(results[0], 3) == 0.532, "AFRQ[0]: got " + str(results[0])
    assert round(results[1], 3) == 0.75, "AFRQ[1]: got " + str(results[1])
    assert round(results[2], 3) == 0.914, "AFRQ[2]: got " + str(results[2])
}

PROB — Introduction to Random Strings

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Computed in logarithms because the raw probability of any particular long string underflows to zero — the same reason the HMM problems decode in log space rather than multiplying probabilities directly.

# Rosalind: PROB — Introduction to Random Strings
# https://rosalind.info/problems/prob/
#
# Given: A DNA string s and an array A of GC contents.
# Return: For each GC content, the common logarithm of the probability that a
# random string with that GC content matches s exactly.

let s = "ACGATACAA"
let gc_contents = [0.129, 0.287, 0.423, 0.476, 0.641, 0.742, 0.783]

# At GC content x, each of G and C appears with probability x/2 and each of
# A and T with probability (1-x)/2. Sum the logs rather than multiplying the
# probabilities, which would underflow on a long string.
let results = gc_contents |> map(|x| {
    let gc_probability = x / 2.0
    let at_probability = (1.0 - x) / 2.0
    range(0, len(s))
        |> map(|i| {
            let base = substr(s, i, 1)
            let is_gc = base == "G" or base == "C"
            log10(if is_gc then gc_probability else at_probability)
        })
        |> sum()
})

let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: -5.737 -5.217 -5.263 -5.36 -5.958 -6.628 -7.009")

fn test_prob_random_string_logs() {
    assert round(results[0], 3) == -5.737, "PROB[0]: got " + str(results[0])
    assert round(results[3], 3) == -5.36, "PROB[3]: got " + str(results[3])
    assert round(results[6], 3) == -7.009, "PROB[6]: got " + str(results[6])
}

LEXF — Enumerating k-mers Lexicographically

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Enumerating every k-mer over an alphabet, which is 4 to the k for DNA. Fine at k=3 and hopeless by k=20, which is why the k-mer problems index sequences rather than enumerate possibilities.

# Rosalind: LEXF — Enumerating k-mers Lexicographically
# https://rosalind.info/problems/lexf/
#
# Given: A collection of at most 10 symbols and a positive integer n <= 10.
# Return: All strings of length n that can be formed from the alphabet,
# ordered lexicographically.

let alphabet = ["A", "C", "G", "T"]
let n = 2

# Grow the set one position at a time. Because the alphabet is already in
# order and each existing prefix is extended in order, the result comes out
# lexicographically sorted without a final sort.
fn words_of_length(symbols, length) {
    range(0, length) |> reduce(
        |acc, _| acc |> flat_map(|prefix| symbols |> map(|s| prefix ++ s)),
        [""]
    )
}

let words = words_of_length(alphabet, n)

println("Result:   " + str(len(words)) + " words, first four: " + (take(words, 4) |> join(" ")))
println("Expected: 16 words, first four: AA AC AG AT")

fn test_lexf_lexicographic_kmers() {
    assert len(words) == 16, "LEXF: got " + str(len(words)) + " words"
    assert words[0] == "AA", "LEXF: first is " + words[0]
    assert words[15] == "TT", "LEXF: last is " + words[15]
    # Bound first: `x |> join(" ") == "..."` binds the comparison to join's
    # second argument rather than to its result.
    let first_four = take(words, 4) |> join(" ")
    assert first_four == "AA AC AG AT", "LEXF: order was " + first_four
}

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"
}

TRAN — Transitions and Transversions

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A random mutation process would give a ratio near 0.5, since there are twice as many ways to make a transversion. Real human data sits near 2.1 genome-wide and near 3 in exomes, largely because methylated cytosine deaminates to thymine. The ratio is therefore a standard quality check: a callset far from those values is reporting sequencing error rather than biology.

# Rosalind: TRAN — Transitions and Transversions
# https://rosalind.info/problems/tran/
#
# Given: Two DNA strings of equal length.
# Return: The transition/transversion ratio.

let s1 = "GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGAAGTACGGGCATCAACCCAGTT"
let s2 = "TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGCGGTACGAGTGTTCCTTTGGGT"

# A transition swaps purine for purine (A<->G) or pyrimidine for pyrimidine
# (C<->T); anything else is a transversion.
fn is_purine(base) { base == "A" or base == "G" }

fn is_transition(a, b) { is_purine(a) == is_purine(b) }

let differing = range(0, len(s1)) |> filter(|i| substr(s1, i, 1) != substr(s2, i, 1))
let transitions = differing |> count_if(|i| is_transition(substr(s1, i, 1), substr(s2, i, 1)))
let transversions = len(differing) - transitions

let ratio = float(transitions) / float(transversions)

println("Result:   " + str(ratio))
println("Expected: 1.21428571429")

fn test_tran_ratio() {
    assert round(ratio, 8) == 1.21428571, "TRAN: got " + str(ratio)
}

SPLC — RNA Splicing

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

What RNA glossed over. A eukaryotic gene is not a contiguous coding sequence — introns are cut out before translation, so the protein comes from the exons joined together. Removing the introns first is the difference between the right protein and nonsense.

# Rosalind: SPLC — RNA Splicing
# https://rosalind.info/problems/splc/
#
# Given: A DNA string s and a collection of introns, all in FASTA format.
# Return: The protein string translated from the exons of s.

let pre_mrna = "ATGGTCTACATAGCTGACAAACAGCACGTAGCAATCGGTCGAATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCGCCTAG"
let introns = [
    "ATCGGTCGAA",
    "ATCGGTCGAGCGTGT"
]

# Remove each intron once, in the order given, then translate what is left.
let coding = introns |> reduce(|seq, intron| replace(seq, intron, ""), pre_mrna)
let peptide = translate(dna(coding))

println("Result:   " ++ str(peptide))
println("Expected: MVYIADKQHVASREAYGHMFKVCA")

fn test_splc_rna_splicing() {
    assert str(peptide) == "MVYIADKQHVASREAYGHMFKVCA", "SPLC: got " ++ str(peptide)
}

CONS — Consensus and Profile

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A profile is the same object BA2's motif search builds and scores against: counts per base per position. The consensus reads off the commonest base at each, which is a useful summary and a lossy one — it discards how strong the preference was.

# Rosalind: CONS — Consensus and Profile
# https://rosalind.info/problems/cons/
#
# Given: A collection of DNA strings of equal length in FASTA format.
# Return: A consensus string, and the profile matrix of base counts per column.

let strings = [
    "ATCCAGCT", "GGGCAACT", "ATGGATCT", "AAGCAACC",
    "TTGGAACT", "ATGCCATT", "ATGGCACT"
]

let bases = ["A", "C", "G", "T"]
# Held as plain strings: str() on a DNA value renders as "DNA(...)", so substr
# would slice the wrapper rather than the sequence.
let width = len(strings[0])

# profile[base][column] — how often each base appears at each position.
let profile = bases |> map(|b| {
    range(0, width) |> map(|i| strings |> count_if(|s| substr(s, i, 1) == b))
})

# The consensus takes whichever base is commonest in each column.
let consensus_string = range(0, width)
    |> map(|i| {
        let column = range(0, 4) |> map(|b| { base: bases[b], n: profile[b][i] })
        let best = column |> sort_by(|e| e.n) |> reverse() |> first()
        best.base
    })
    |> join("")

println("Result:   " + consensus_string)
println("Expected: ATGCAACT")
range(0, 4) |> each(|b| println("  " + bases[b] + ": " + (profile[b] |> map(|n| str(n)) |> join(" "))))

fn test_cons_consensus_and_profile() {
    assert consensus_string == "ATGCAACT", "CONS: got " + consensus_string
    assert profile[0][0] == 5, "CONS: A count at column 1 was " + str(profile[0][0])
    assert profile[3][0] == 1, "CONS: T count at column 1 was " + str(profile[3][0])
}

GRPH — Overlap Graphs

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Overlap graphs make reads the nodes, which makes assembly a Hamiltonian path — NP-hard. DBRU makes reads the edges instead and gets an Eulerian path, solvable in linear time. Same data, different graph, and the whole reason modern assemblers use de Bruijn graphs.

# Rosalind: GRPH — Overlap Graphs
# https://rosalind.info/problems/grph/
#
# Given: A collection of DNA strings in FASTA format.
# Return: The adjacency list of the overlap graph O(3), where an edge s -> t
# means the last 3 characters of s equal the first 3 of t, and s != t.

let k = 3
let records = [
    { id: "Rosalind_0498", seq: "AAATAAA" },
    { id: "Rosalind_2391", seq: "AAATTTT" },
    { id: "Rosalind_2323", seq: "TTTTCCC" },
    { id: "Rosalind_0442", seq: "AAATCCC" },
    { id: "Rosalind_5013", seq: "GGGTGGG" }
]

fn suffix_of(s, width) { substr(s, len(s) - width, width) }
fn prefix_of(s, width) { substr(s, 0, width) }

let edge_list = records |> flat_map(|a| {
    records
        |> filter(|b| a.id != b.id and suffix_of(a.seq, k) == prefix_of(b.seq, k))
        |> map(|b| a.id ++ " " ++ b.id)
})

println("Result:")
edge_list |> each(|e| println("  " + e))
println("Expected: 3 edge_list — 0498->2391, 0498->0442, 2391->2323")

fn test_grph_overlap_graph() {
    assert len(edge_list) == 3, "GRPH: got " + str(len(edge_list)) + " edge_list"
    assert edge_list |> contains("Rosalind_0498 Rosalind_2391"), "GRPH: missing 0498->2391"
    assert edge_list |> contains("Rosalind_0498 Rosalind_0442"), "GRPH: missing 0498->0442"
    assert edge_list |> contains("Rosalind_2391 Rosalind_2323"), "GRPH: missing 2391->2323"
}

REVP — Locating Restriction Sites

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Restriction enzymes cut at reverse palindromes because they bind as symmetric dimers, one subunit per strand, so the site reads the same on both. That is why these particular sequences matter and why they are always even in length.

# Rosalind: REVP — Locating Restriction Sites
# https://rosalind.info/problems/revp/
#
# Given: A DNA string of length at most 1 kbp.
# Return: The position and length of every reverse palindrome of length
# between 4 and 12, as 1-based positions.

let s = "TCAATGCATGCGGGTCTATATGCAT"

# A reverse palindrome reads the same as its own reverse complement, so a
# window qualifies when it equals reverse_complement(window).
fn is_reverse_palindrome(window) {
    str(reverse_complement(dna(window))) == window
}

let found = range(4, 13) |> flat_map(|width| {
    range(0, len(s) - width + 1)
        |> filter(|i| is_reverse_palindrome(substr(s, i, width)))
        |> map(|i| { position: i + 1, length: width })
}) |> sort_by(|hit| hit.position)

println("Result:   " + str(len(found)) + " palindromes")
found |> each(|hit| println("  " + str(hit.position) + " " + str(hit.length)))
println("Expected: 8 — (4,6) (5,4) (6,6) (7,4) (17,4) (18,4) (20,6) (21,4)")

fn test_revp_reverse_palindromes() {
    assert len(found) == 8, "REVP: got " + str(len(found))
    assert found[0].position == 4 and found[0].length == 6, "REVP: first was " + str(found[0])
}

TREE — Completing a Tree

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A tree on n nodes has exactly n-1 edges, so the answer is that count minus the edges already present, with no search involved. The point is recognising the invariant rather than exploring the graph.

# Rosalind: TREE — Completing a Tree
# https://rosalind.info/problems/tree/
#
# Given: A positive integer n and an adjacency list of a graph on n nodes that
# forms a forest.
# Return: The minimum number of edges needed to produce a tree.

let n = 10
let edge_list = [[1,2],[2,8],[4,10],[5,9],[6,10],[7,9]]

# A tree on n nodes has exactly n-1 edges, and a forest with e edges has
# n-e components — so joining them needs (n-1)-e more.
let result = (n - 1) - len(edge_list)

println("Result:   " + str(result))
println("Expected: 3")

fn test_tree_edges_to_add() {
    assert result == 3, "TREE: got " + str(result)
}

INOD — Counting Phylogenetic Ancestors

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

An unrooted binary tree with n leaves always has n-2 internal nodes, whatever its shape. That fixed relationship is what lets EUBT count trees by insertion and CNTQ count quartets without ever examining a topology.

# Rosalind: INOD — Counting Phylogenetic Ancestors
# https://rosalind.info/problems/inod/
#
# Given: A positive integer n (3 <= n <= 10000).
# Return: The number of internal nodes of any unrooted binary tree with n leaves.

let n = 4

# Every internal node of an unrooted binary tree has degree 3. Counting edge
# endpoints two ways gives internal = n - 2, independent of the shape.
let result = n - 2

println("Result:   " + str(result))
println("Expected: 2")

fn test_inod_internal_nodes() {
    assert result == 2, "INOD: got " + str(result)
}

SIGN — Enumerating Oriented Gene Orderings

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Signed permutations: 2^n times n! of them, because each gene may also be flipped. The sign is not decoration — a reversed gene reads on the opposite strand, which is exactly what REAR and SORT measure.

# Rosalind: SIGN — Enumerating Oriented Gene Orderings
# https://rosalind.info/problems/sign/
#
# Given: A positive integer n <= 6.
# Return: The total number of signed permutations of length n, followed by a
# list of all such permutations.

let n = 2

fn permutations_of(items) {
    if len(items) <= 1 {
        [items]
    } else {
        range(0, len(items)) |> flat_map(|i| {
            let rest = concat(take(items, i), drop(items, i + 1))
            permutations_of(rest) |> map(|p| concat([items[i]], p))
        })
    }
}

# Every permutation admits 2^n sign assignments, applied independently.
fn sign_patterns(width) {
    range(0, width) |> reduce(
        |acc, _| acc |> flat_map(|p| [concat(p, [1]), concat(p, [-1])]),
        [[]]
    )
}

let signed = permutations_of(range(1, n + 1)) |> flat_map(|p| {
    sign_patterns(n) |> map(|signs| range(0, n) |> map(|i| p[i] * signs[i]))
})

println("Result:   " + str(len(signed)))
println("Expected: 8")
signed |> each(|p| println("  " + (p |> map(|x| str(x)) |> join(" "))))

fn test_sign_signed_permutations() {
    assert len(signed) == 8, "SIGN: got " + str(len(signed))
    let distinct = signed |> map(|p| p |> map(|x| str(x)) |> join(" ")) |> unique()
    assert len(distinct) == 8, "SIGN: duplicates present"
}

LGIS — Longest Increasing Subsequence

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Longest increasing subsequence, which in a genomic setting is how conserved order is recovered from a comparison of two genomes: positions that stay in ascending order form a synteny block, and the rest is rearrangement.

# Rosalind: LGIS — Longest Increasing Subsequence
# https://rosalind.info/problems/lgis/
#
# Given: A positive integer n and a permutation of length n.
# Return: A longest increasing subsequence, then a longest decreasing one.

let permutation = [5, 1, 4, 2, 3]

# O(n^2) dynamic programme: best[i] is the length of the longest run ending at
# i, and prev[i] remembers which index it came from so the run can be rebuilt.
fn longest_run(values, increasing) {
    let n = len(values)
    let best = range(0, n) |> map(|_| 1)
    let prev = range(0, n) |> map(|_| -1)

    let i = 1
    while i < n {
        let j = 0
        while j < i {
            let ordered = if increasing then values[j] < values[i] else values[j] > values[i]
            if ordered and best[j] + 1 > best[i] {
                best = set_at(best, i, best[j] + 1)
                prev = set_at(prev, i, j)
            }
            j = j + 1
        }
        i = i + 1
    }

    let last = 0
    let k = 1
    while k < n {
        if best[k] > best[last] then last = k
        k = k + 1
    }

    let run = []
    let at = last
    while at >= 0 {
        run = concat([values[at]], run)
        at = prev[at]
    }
    run
}

fn set_at(list, index, value) {
    range(0, len(list)) |> map(|i| if i == index then value else list[i])
}

let increasing = longest_run(permutation, true)
let decreasing = longest_run(permutation, false)

println("Increasing: " + (increasing |> map(|x| str(x)) |> join(" ")))
println("Decreasing: " + (decreasing |> map(|x| str(x)) |> join(" ")))
println("Expected:   1 2 3  /  5 4 2")

fn test_lgis_longest_runs() {
    assert len(increasing) == 3, "LGIS: increasing length " + str(len(increasing))
    assert len(decreasing) == 3, "LGIS: decreasing length " + str(len(decreasing))
}

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 + "'"
}

PDST — Creating a Distance Matrix

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Pairwise p-distance: the fraction of positions at which two sequences differ. It is the crude measure HAMM warns about, applied to every pair, and it feeds the tree-building problems that follow.

# Rosalind: PDST — Creating a Distance Matrix
# https://rosalind.info/problems/pdst/
#
# Given: A collection of DNA strings of equal length in FASTA format.
# Return: The matrix D of p-distances, where D[i][j] is the proportion of
# positions at which strings i and j differ.

let strings = [
    "TTTCCATTTA",
    "GATTCATTTC",
    "TTTCCATTTT",
    "GTTCCATTTA"
]

fn p_distance(a, b) {
    let differing = range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1))
    float(differing) / float(len(a))
}

let distance_grid = strings |> map(|a| strings |> map(|b| p_distance(a, b)))

println("Result:")
distance_grid |> each(|row| println("  " + (row |> map(|d| str(round(d, 5))) |> join(" "))))
println("Expected first row: 0.00000 0.40000 0.10000 0.10000")

fn test_pdst_distance_matrix() {
    assert round(distance_grid[0][0], 5) == 0.0, "PDST: diagonal not zero"
    assert round(distance_grid[0][1], 5) == 0.4, "PDST: [0][1] was " + str(distance_grid[0][1])
    assert round(distance_grid[0][2], 5) == 0.1, "PDST: [0][2] was " + str(distance_grid[0][2])
    assert round(distance_grid[0][3], 5) == 0.1, "PDST: [0][3] was " + str(distance_grid[0][3])
    assert round(distance_grid[1][0], 5) == round(distance_grid[0][1], 5), "PDST: not symmetric"
}

ASMQ — Assessing Assembly Quality with N50 and N75

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

N50 is the contig length at which half the assembly sits in contigs that long or longer. It is the standard summary and a gameable one: joining contigs wrongly raises it, so a high N50 is evidence of a long assembly rather than a correct one.

# Rosalind: ASMQ — Assessing Assembly Quality with N50 and N75
# https://rosalind.info/problems/asmq/
#
# Given: A collection of at most 1000 DNA strings.
# Return: N50 and N75 for the collection.

let contigs = [
    "GATTACA",
    "TACTACTAC",
    "ATTGAT",
    "GAAGA"
]

# NXX is the length of the shortest contig in the set of longest contigs that
# together cover at least XX% of the assembly.
fn n_statistic(lengths, percent) {
    let sorted = lengths |> sort() |> reverse()
    let target = float(sum(sorted)) * percent / 100.0
    let running = 0.0
    let answer = 0
    let i = 0
    while i < len(sorted) and answer == 0 {
        running = running + float(sorted[i])
        if running >= target then answer = sorted[i]
        i = i + 1
    }
    answer
}

let lengths = contigs |> map(|c| len(c))
let n50_value = n_statistic(lengths, 50.0)
let n75 = n_statistic(lengths, 75.0)

println("Result:   " + str(n50_value) + " " + str(n75))
println("Expected: 7 6")

fn test_asmq_n50_and_n75() {
    assert n50_value == 7, "ASMQ: N50 was " + str(n50_value)
    assert n75 == 6, "ASMQ: N75 was " + str(n75)
}

ORF — Open Reading Frames

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Six frames, not three — the gene may sit on either strand, and there is no way to tell from the sequence which. Each frame is scanned from every start codon to the first stop.

# Rosalind: ORF — Open Reading Frames
# https://rosalind.info/problems/orf/
#
# Given: A DNA string s of length at most 1 kbp.
# Return: Every distinct protein that can be translated from an ORF of s,
# considering both strands.

let s = "AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG"

let stop_codons = ["TAA", "TAG", "TGA"]

# Find the in-frame stop, then translate the segment in one call. An ORF that
# runs off the end without a stop does not count.
fn orf_from(strand, start) {
    let i = start
    let stop_at = -1
    while i + 3 <= len(strand) and stop_at < 0 {
        if stop_codons |> contains(substr(strand, i, 3)) then stop_at = i
        i = i + 3
    }
    if stop_at < 0 {
        ""
    } else {
        str(translate(dna(substr(strand, start, stop_at - start))))
    }
}

fn orfs_of(strand) {
    range(0, len(strand) - 2)
        |> filter(|i| substr(strand, i, 3) == "ATG")
        |> map(|i| orf_from(strand, i))
        |> filter(|p| p != "")
}

let reverse_strand = str(reverse_complement(dna(s)))

let proteins = concat(orfs_of(s), orfs_of(reverse_strand)) |> unique()

println("Result:   " + str(len(proteins)) + " distinct proteins")
proteins |> each(|p| println("  " + p))
println("Expected: 4 — MLLGSFRLIPKETLIQVAGSSPCNLS, M, MGMTPRLGLESLLE, MTPRLGLESLLE")

fn test_orf_open_reading_frames() {
    assert len(proteins) == 4, "ORF: got " + str(len(proteins))
    assert proteins |> contains("MLLGSFRLIPKETLIQVAGSSPCNLS"), "ORF: missing the long reverse-strand protein"
    assert proteins |> contains("MGMTPRLGLESLLE"), "ORF: missing MGMTPRLGLESLLE"
    assert proteins |> contains("MTPRLGLESLLE"), "ORF: missing MTPRLGLESLLE"
    assert proteins |> contains("M"), "ORF: missing the single-residue ORF"
}

LEXV — Ordering Strings of Varying Length Lexicographically

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Ordering strings of unequal length, where a prefix sorts before anything extending it. Ordinary lexicographic order on padded strings gets this wrong, which is the whole exercise.

# Rosalind: LEXV — Ordering Strings of Varying Length Lexicographically
# https://rosalind.info/problems/lexv/
#
# Given: An ordered alphabet of at most 10 symbols and a positive integer n.
# Return: All strings of length at most n from that alphabet, ordered by the
# alphabet's own order rather than by ASCII.

let alphabet = ["D", "N", "A"]
let n = 3

# Pre-order walk of the trie: emit a string, then everything that extends it.
# That places "D" before "DD", which is what "varying length" ordering means
# here — a prefix sorts before anything built on it.
fn extensions(prefix, symbols, remaining) {
    if remaining == 0 {
        []
    } else {
        symbols |> flat_map(|s| {
            let word = prefix ++ s
            concat([word], extensions(word, symbols, remaining - 1))
        })
    }
}

let words = extensions("", alphabet, n)

println("Result:   " + str(len(words)) + " strings")
println("First 8:  " + (take(words, 8) |> join(" ")))
println("Expected: 39 strings, starting D DD DDD DDN DDA DN DND DNN")

fn test_lexv_varying_length_order() {
    assert len(words) == 39, "LEXV: got " + str(len(words))
    let first_eight = take(words, 8) |> join(" ")
    assert first_eight == "D DD DDD DDN DDA DN DND DNN", "LEXV: order was " + first_eight
    assert words[38] == "AAA", "LEXV: last was " + words[38]
}

SETO — Introduction to Set Operations

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Set operations as groundwork. Worth having explicitly because several later problems phrase their answer as a complement or an intersection, and getting the universe wrong is a common way to be quietly off by a few elements.

# Rosalind: SETO — Introduction to Set Operations
# https://rosalind.info/problems/seto/
#
# Given: A positive integer n and two subsets A and B of {1, ..., n}.
# Return: Their union, intersection, both differences, and both complements.

let n = 10
let a = [1, 2, 3, 4, 5]
let b = [2, 8, 5, 10]

let universe = range(1, n + 1)

fn union_of(x, y) { concat(x, y) |> unique() |> sort() }
fn intersection_of(x, y) { x |> filter(|e| y |> contains(e)) |> sort() }
fn difference_of(x, y) { x |> filter(|e| !(y |> contains(e))) |> sort() }

let combined_items = union_of(a, b)
let shared_items = intersection_of(a, b)
let a_minus_b = difference_of(a, b)
let b_minus_a = difference_of(b, a)
let complement_a = difference_of(universe, a)
let complement_b = difference_of(universe, b)

fn show(label, s) { println("  " + label + ": " + (s |> map(|e| str(e)) |> join(" "))) }

println("Result:")
show("A combined_items B     ", combined_items)
show("A inter B     ", shared_items)
show("A - B         ", a_minus_b)
show("B - A         ", b_minus_a)
show("complement A  ", complement_a)
show("complement B  ", complement_b)

fn joined(s) { s |> map(|e| str(e)) |> join(" ") }

fn test_seto_set_operations() {
    assert joined(combined_items) == "1 2 3 4 5 8 10", "SETO combined_items: " + joined(combined_items)
    assert joined(shared_items) == "2 5", "SETO shared_items: " + joined(shared_items)
    assert joined(a_minus_b) == "1 3 4", "SETO A-B: " + joined(a_minus_b)
    assert joined(b_minus_a) == "8 10", "SETO B-A: " + joined(b_minus_a)
    assert joined(complement_a) == "6 7 8 9 10", "SETO ~A: " + joined(complement_a)
    assert joined(complement_b) == "1 3 4 6 7 9", "SETO ~B: " + joined(complement_b)
}

LCSQ — Finding a Shared Spliced Motif

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Any longest common subsequence is a valid answer, so the assertion checks the length and that the motif really is a subsequence of both strings.

# Rosalind: LCSQ — Finding a Shared Spliced Motif
# https://rosalind.info/problems/lcsq/
#
# Given: Two DNA strings s and t.
# Return: A longest common subsequence of s and t.

let s = "AACCTTGG"
let t = "ACACTGTGA"

# Standard LCS table, then walk it backwards to rebuild one optimal answer.
# Any longest common subsequence is accepted, so the traceback's tie-breaking
# does not matter.
fn lcs_of(a, b) {
    let rows = len(a) + 1
    let cols = len(b) + 1
    let table = range(0, rows) |> map(|_| range(0, cols) |> map(|_| 0))

    let i = 1
    while i < rows {
        let ai = substr(a, i - 1, 1)
        let row = [0]
        let j = 1
        while j < cols {
            let value = if ai == substr(b, j - 1, 1) {
                table[i - 1][j - 1] + 1
            } else {
                max([table[i - 1][j], row[j - 1]])
            }
            row = push(row, value)
            j = j + 1
        }
        table = set_row(table, i, row)
        i = i + 1
    }

    # Traceback from the bottom-right corner.
    let result = ""
    let x = len(a)
    let y = len(b)
    while x > 0 and y > 0 {
        if substr(a, x - 1, 1) == substr(b, y - 1, 1) {
            result = substr(a, x - 1, 1) ++ result
            x = x - 1
            y = y - 1
        } else {
            if table[x - 1][y] >= table[x][y - 1] then x = x - 1 else y = y - 1
        }
    }
    result
}

fn set_row(table, index, row) {
    range(0, len(table)) |> map(|i| if i == index then row else table[i])
}

fn is_subsequence(needle, haystack) {
    let at = 0
    let i = 0
    while i < len(haystack) and at < len(needle) {
        if substr(haystack, i, 1) == substr(needle, at, 1) then at = at + 1
        i = i + 1
    }
    at == len(needle)
}

let motif = lcs_of(s, t)

println("Result:   " + motif + " (length " + str(len(motif)) + ")")
println("Expected: a length-6 common subsequence, such as AACTGG")

fn test_lcsq_shared_spliced_motif() {
    assert len(motif) == 6, "LCSQ: length " + str(len(motif)) + " for '" + motif + "'"
    assert is_subsequence(motif, s), "LCSQ: '" + motif + "' is not a subsequence of s"
    assert is_subsequence(motif, t), "LCSQ: '" + motif + "' is not a subsequence of t"
}

KMP — Speeding Up Motif Finding

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The failure array records, for each prefix, the longest proper prefix that is also a suffix — so a mismatch can resume without re-reading the text. That is what makes matching linear rather than quadratic, and it is the same insight the trie and suffix-array problems generalise.

# Rosalind: KMP — Speeding Up Motif Finding
# https://rosalind.info/problems/kmp/
#
# Given: A DNA string s.
# Return: The failure array of s: for each prefix, the length of the longest
# proper prefix that is also a suffix of it.

let s = "CAGCATGGTATCACAGCAGAG"

# The Knuth-Morris-Pratt construction: extend the previous border where the
# next character agrees, otherwise fall back through shorter borders.
let failure = [0]
let border = 0
let i = 1
while i < len(s) {
    while border > 0 and substr(s, i, 1) != substr(s, border, 1) {
        border = failure[border - 1]
    }
    if substr(s, i, 1) == substr(s, border, 1) then border = border + 1
    failure = push(failure, border)
    i = i + 1
}

let formatted = failure |> map(|v| str(v)) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0 0 0 1 2 0 0 0 0 0 0 1 2 1 2 3 4 5 3 0 0")

fn test_kmp_failure_array() {
    assert formatted == "0 0 0 1 2 0 0 0 0 0 0 1 2 1 2 3 4 5 3 0 0", "KMP: got " + formatted
}

DBRU — Constructing a De Bruijn Graph

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Reads become edges rather than nodes, which turns assembly from a Hamiltonian path problem into an Eulerian one. GRPH shows the alternative; this is the construction real assemblers are built on.

# Rosalind: DBRU — Constructing a De Bruijn Graph
# https://rosalind.info/problems/dbru/
#
# Given: A collection of up to 1000 DNA strings of equal length, all distinct.
# Return: The adjacency list of the De Bruijn graph on S union S_rc, where each
# k-mer contributes an edge from its prefix to its suffix.

let reads = ["TGAT", "CATG", "TCAT", "ATGC", "CATC", "CATC"]

# The graph is built on the reads together with their reverse complements,
# as a set — the duplicate CATC and the shared CATG collapse.
let with_reverse = reads |> flat_map(|r| [r, str(reverse_complement(dna(r)))]) |> unique()

let edge_list = with_reverse
    |> map(|kmer| substr(kmer, 0, len(kmer) - 1) ++ " -> " ++ substr(kmer, 1, len(kmer) - 1))
    |> unique()
    |> sort()

println("Result:   " + str(len(edge_list)) + " edge_list")
edge_list |> each(|e| println("  (" + replace(e, " -> ", ", ") + ")"))
println("Expected: 9 edge_list, starting (ATC, TCA) (ATG, TGA) (ATG, TGC)")

fn test_dbru_de_bruijn_graph() {
    assert len(edge_list) == 9, "DBRU: got " + str(len(edge_list))
    assert edge_list[0] == "ATC -> TCA", "DBRU: first edge " + edge_list[0]
    assert edge_list |> contains("GCA -> CAT"), "DBRU: missing GCA -> CAT"
    assert edge_list |> contains("TGA -> GAT"), "DBRU: missing TGA -> GAT"
}

ASPC — Introduction to Alternative Splicing

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Summing binomial coefficients, motivated by alternative splicing: one gene yields many proteins depending on which exons are kept. The count grows fast enough that the modulus is doing real work.

# Rosalind: ASPC — Introduction to Alternative Splicing
# https://rosalind.info/problems/aspc/
#
# Given: Positive integers n and m with n >= m.
# Return: The sum of C(n, k) for all k from m to n, modulo 1,000,000.

let n = 6
let m = 3
let modulus = 1000000

# Build Pascal's triangle row by row, reducing as we go, so nothing overflows
# even when n reaches 2000.
fn binomials(width, md) {
    range(0, width) |> reduce(
        |row_values, _| {
            let extended = concat([0], row_values)
            range(0, len(row_values) + 1) |> map(|i| {
                let left = extended[i]
                let right = if i < len(row_values) then row_values[i] else 0
                (left + right) % md
            })
        },
        [1]
    )
}

let row_values = binomials(n, modulus)
# Bound before the modulo: `|> sum() % modulus` reads the modulo as another
# argument to sum().
let total = range(m, n + 1) |> map(|k| row_values[k]) |> sum()
let result = total % modulus

println("Result:   " + str(result))
println("Expected: 42")

fn test_aspc_alternative_splicing() {
    assert result == 42, "ASPC: got " + str(result)
}

CORR — Error Correction in Reads

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Assumes a correct read appears at least twice and an erroneous one exactly once — which is what read depth buys and why coverage matters. Reverse complements count as the same read, since strand is not recorded.

# Rosalind: CORR — Error Correction in Reads
# https://rosalind.info/problems/corr/
#
# Given: A collection of reads of equal length. Each read either appears in the
# dataset at least twice (counting reverse complements), or is a single-symbol
# error away from exactly one such read.
# Return: Every correction, as "wrong->right".

let reads = ["TCATC", "TTCAT", "TCATC", "TGAAA", "GAGGA", "TTTCA", "ATCAA", "TTGAT", "TTTCC"]

fn revcomp(read) { str(reverse_complement(dna(read))) }

fn hamming_of(a, b) {
    range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1))
}

# A read is correct when it, plus its reverse complement, appears at least
# twice across the dataset.
fn occurrences(read) {
    let direct = reads |> count_if(|r| r == read)
    let reversed = reads |> count_if(|r| r == revcomp(read))
    direct + reversed
}

let correct = reads |> filter(|r| occurrences(r) >= 2) |> unique()
let wrong = reads |> filter(|r| occurrences(r) < 2) |> unique()

# Each incorrect read is one substitution from exactly one correct read, in
# either orientation.
let corrections = wrong |> map(|bad| {
    let targets = correct |> flat_map(|good| [good, revcomp(good)])
    let repaired = targets |> filter(|t| hamming_of(bad, t) == 1) |> first()
    bad ++ "->" ++ repaired
})

println("Result:")
corrections |> each(|c| println("  " + c))
println("Expected: TTCAT->TTGAT, GAGGA->GATGA, TTTCC->TTTCA")

fn test_corr_error_correction() {
    assert len(corrections) == 3, "CORR: got " + str(len(corrections))
    assert corrections |> contains("TTCAT->TTGAT"), "CORR: missing TTCAT->TTGAT"
    assert corrections |> contains("GAGGA->GATGA"), "CORR: missing GAGGA->GATGA"
    assert corrections |> contains("TTTCC->TTTCA"), "CORR: missing TTTCC->TTTCA"
}

SCSP — Interleaving Two Motifs

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Any shortest common supersequence is valid, so the assertion checks the length and that both inputs are subsequences of the result.

# Rosalind: SCSP — Interleaving Two Motifs
# https://rosalind.info/problems/scsp/
#
# Given: Two DNA strings s and t.
# Return: A shortest common supersequence of s and t.

let s = "ATCTGAT"
let t = "TGCATA"

fn set_row(grid, index, row) {
    range(0, len(grid)) |> map(|i| if i == index then row else grid[i])
}

# The shortest common supersequence is built from the longest common
# subsequence: walk both strings, emitting shared characters once and the rest
# as they come.
fn lcs_table(a, b) {
    let grid = range(0, len(a) + 1) |> map(|_| range(0, len(b) + 1) |> map(|_| 0))
    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        let row = [0]
        let j = 1
        while j <= len(b) {
            let value = if ai == substr(b, j - 1, 1) {
                grid[i - 1][j - 1] + 1
            } else {
                max([grid[i - 1][j], row[j - 1]])
            }
            row = push(row, value)
            j = j + 1
        }
        grid = set_row(grid, i, row)
        i = i + 1
    }
    grid
}

let grid = lcs_table(s, t)

let result = ""
let x = len(s)
let y = len(t)
while x > 0 and y > 0 {
    if substr(s, x - 1, 1) == substr(t, y - 1, 1) {
        result = substr(s, x - 1, 1) ++ result
        x = x - 1
        y = y - 1
    } else {
        if grid[x - 1][y] >= grid[x][y - 1] {
            result = substr(s, x - 1, 1) ++ result
            x = x - 1
        } else {
            result = substr(t, y - 1, 1) ++ result
            y = y - 1
        }
    }
}
result = substr(s, 0, x) ++ substr(t, 0, y) ++ result

fn is_subsequence(needle, haystack) {
    let at = 0
    let i = 0
    while i < len(haystack) and at < len(needle) {
        if substr(haystack, i, 1) == substr(needle, at, 1) then at = at + 1
        i = i + 1
    }
    at == len(needle)
}

println("Result:   " + result + " (length " + str(len(result)) + ")")
println("Expected: a length-9 supersequence, such as ATGCATGAT")

fn test_scsp_shortest_common_supersequence() {
    assert len(result) == 9, "SCSP: length " + str(len(result)) + " for '" + result + "'"
    assert is_subsequence(s, result), "SCSP: s is not a subsequence of the result"
    assert is_subsequence(t, result), "SCSP: t is not a subsequence of the result"
}

EVAL — Expected Number of Restriction Sites

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The expected number of times a motif appears in a random sequence, which is what makes an observed count interesting or unremarkable. Without a baseline, finding a site says nothing at all.

# Rosalind: EVAL — Expected Number of Restriction Sites
# https://rosalind.info/problems/eval/
#
# Given: A positive integer n, a DNA string s, and an array A of GC contents.
# Return: For each GC content, the expected number of times s appears as a
# substring of a random string of length n.

let n = 10
let s = "AG"
let gc_contents = [0.25, 0.5, 0.75]

# Each of the (n - |s| + 1) starting positions carries the same probability of
# spelling s, and expectation is additive whether or not those events overlap.
let positions = n - len(s) + 1

let results = gc_contents |> map(|x| {
    # No product() builtin, so fold the multiplication.
    let probability = range(0, len(s)) |> reduce(|acc, i| {
        let base = substr(s, i, 1)
        let is_gc = base == "G" or base == "C"
        acc * (if is_gc then x / 2.0 else (1.0 - x) / 2.0)
    }, 1.0)
    float(positions) * probability
})

let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0.422 0.563 0.422")

fn test_eval_expected_sites() {
    assert round(results[0], 3) == 0.422, "EVAL[0]: got " + str(results[0])
    assert round(results[1], 3) == 0.563, "EVAL[1]: got " + str(results[1])
    assert round(results[2], 3) == 0.422, "EVAL[2]: got " + str(results[2])
}

LONG — Genome Assembly as Shortest Superstring

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Shortest superstring is NP-hard in general. This is solvable only because the problem guarantees every pair overlaps by more than half their length, which makes the correct overlap unique and a greedy merge safe. GREP shows what happens when that guarantee is dropped.

# Rosalind: LONG — Genome Assembly as Shortest Superstring
# https://rosalind.info/problems/long/
#
# Given: At most 50 DNA strings of equal length, where every pair overlaps by
# more than half their length.
# Return: The shortest superstring containing all of them.

let reads = ["ATTAGACCTG", "CCTGCCGGAA", "AGACCTGCCG", "GCCGGAATAC"]

# Longest suffix of a that is also a prefix of b, considering only overlaps
# longer than half — which the problem guarantees is unique.
fn overlap_length(a, b) {
    let limit = min([len(a), len(b)])
    let best = 0
    let width = limit
    while width > limit / 2 and best == 0 {
        if substr(a, len(a) - width, width) == substr(b, 0, width) then best = width
        width = width - 1
    }
    best
}

# Repeatedly glue on whichever remaining read overlaps the growing assembly.
let remaining = drop(reads, 1)
let assembled = reads[0]

while len(remaining) > 0 {
    let joined = false
    let next_remaining = []
    for candidate in remaining {
        if !joined and overlap_length(assembled, candidate) > 0 {
            assembled = assembled ++ substr(candidate, overlap_length(assembled, candidate), len(candidate) - overlap_length(assembled, candidate))
            joined = true
        } else {
            if !joined and overlap_length(candidate, assembled) > 0 {
                assembled = substr(candidate, 0, len(candidate) - overlap_length(candidate, assembled)) ++ assembled
                joined = true
            } else {
                next_remaining = push(next_remaining, candidate)
            }
        }
    }
    remaining = next_remaining
}

println("Result:   " + assembled)
println("Expected: ATTAGACCTGCCGGAATAC")

fn test_long_shortest_superstring() {
    assert assembled == "ATTAGACCTGCCGGAATAC", "LONG: got " + assembled
}

KMER — k-Mer Composition

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The official answer is a 256-entry array. Rather than restate it, the assertion checks the ordering endpoints, that the counts total the number of windows, and cross-checks an entry by independent counting — which is what an ordering or lookup mistake would break.

# Rosalind: KMER — k-Mer Composition
# https://rosalind.info/problems/kmer/
#
# Given: A DNA string s.
# Return: The 4-mer composition of s: the count of every one of the 256
# possible 4-mers, in lexicographic order.

let s = "CTTCGAAAGTTTGGGCCGAGTCTTATAGTCGAATAGCCTCCATATTGGCCTCTATCTGAGCTCTACGCACGCCAGCTCCCGCCCTGGCACTGCTTGATATTGGGTCCTATCTGACCGGACTTGACGTACCTGGGCTGATATTGCCGATCTGATTTGTCCAAAATTGACCATAGTTAAATGGCGGAATGGGCTTTGCCTCGCACCGTAGGCGATTCCCAAGTGCAGGCAAGGATATTGCGAAATATCTTCCAGAGGGCGGAAATTGCTATTCATCCTATTTGGCTTGCTCTTAAGATCACCGATCTAGCTTCGAATCGCCATGCACAGTTTGGGCTCATCGGAAGTCTCGGGCTGGTTTAATTGTGATTGGGCCTTAGCTCTCGGCTGATCTTCGTCACGATATCGCTAGATCGGCTTGACTTCGCATTGGTCTGATTTTTAAGATCTTAGGCTTGACGGCGATCTCGCTCGATCCTGCCTGATAGGCTAGCTCGATTCTCATTAAGCTCTAGCCTTCGGCCTCTGGCTCTTTAAGCGATGCTGCACGATATCTATCGGACTGCTCTTGCTCGATATCGACTGCTTGACGGCTTGATATCGCTCTTGGCTTGCTTGACTCGCTCTTGATCGCTCTAGCTAGCTCGCTCTGATCTTGCACTTGCTCGGCTCTTGCTCGGCTCTTAGCTCGCTCGATCGCTCTAGCTCGCTCGATCTTGCTCGATCTTGCACTGCTCGATCTGCACTTGCTCGATCTTGCTCGATCTAGCT"
let k = 4

let alphabet = ["A", "C", "G", "T"]

# Every 4-mer in lexicographic order, so the counts line up with the ordering
# the problem asks for.
let all_kmers = range(0, k) |> reduce(
    |acc, _| acc |> flat_map(|prefix| alphabet |> map(|b| prefix ++ b)),
    [""]
)

# Count the windows of s once, then look each k-mer up.
let window_list = range(0, len(s) - k + 1) |> map(|i| substr(s, i, k))
let counts = all_kmers |> map(|kmer| window_list |> count_if(|w| w == kmer))

println("Result:   " + str(len(counts)) + " counts, first 12: " + (take(counts, 12) |> map(|c| str(c)) |> join(" ")))
println("Total:    " + str(sum(counts)) + " window_list (expected " + str(len(s) - k + 1) + ")")

fn test_kmer_composition() {
    assert len(counts) == 256, "KMER: got " + str(len(counts)) + " entries"
    assert all_kmers[0] == "AAAA", "KMER: order starts at " + all_kmers[0]
    assert all_kmers[255] == "TTTT", "KMER: order ends at " + all_kmers[255]
    # Every window is counted exactly once, so the counts must total the number
    # of windows — the check that catches an ordering or lookup mistake.
    assert sum(counts) == len(s) - k + 1, "KMER: counts total " + str(sum(counts))
    # Cross-check one entry by counting it independently.
    let aaaa = window_list |> count_if(|w| w == "AAAA")
    assert counts[0] == aaaa, "KMER: AAAA count disagrees"
}

EBIN — Wright-Fisher's Expected Behavior

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Expected allele counts under Wright-Fisher, where expectation is linear and so the answer is a sum of binomial means. WFMD gives the distribution this summarises.

# Rosalind: EBIN — Wright-Fisher's Expected Behavior
# https://rosalind.info/problems/ebin/
#
# Given: A positive integer n and an array P of probabilities.
# Return: For each probability, the expected number of successes in n
# independent Bernoulli trials.

let n = 17
let probabilities = [0.1, 0.2, 0.3]

# The expectation of a binomial is n*p, whether or not the trials interact —
# expectation is linear, so no summation over outcomes is needed.
let results = probabilities |> map(|p| float(n) * p)
let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 1.7 3.4 5.1")

fn test_ebin_expected_successes() {
    assert round(results[0], 3) == 1.7, "EBIN[0]: got " + str(results[0])
    assert round(results[1], 3) == 3.4, "EBIN[1]: got " + str(results[1])
    assert round(results[2], 3) == 5.1, "EBIN[2]: got " + str(results[2])
}

SPEC — Inferring Protein from Spectrum

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The answer is derived from the spectrum's consecutive differences (186.079 W, 131.040 M, 128.059 Q, 71.037 A) rather than restated from memory, and the assertion re-checks each residue against its own gap.

# Rosalind: SPEC — Inferring Protein from Spectrum
# https://rosalind.info/problems/spec/
#
# Given: A list L of n masses forming the prefix spectrum of a protein.
# Return: A protein string whose prefix spectrum is L.

let spectrum = [3524.8542, 3710.9335, 3841.974, 3970.0326, 4041.0697]

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}

# Consecutive masses in a prefix spectrum differ by exactly one residue, so
# each gap identifies the residue nearest it. Leucine and isoleucine share a
# mass; the ordering of `keys` decides which name is reported.
fn residue_for(gap) {
    let names = keys(monoisotopic)
    let scored = names |> map(|name| { name: name, error: abs(monoisotopic[name] - gap) })
    let best = scored |> sort_by(|e| e.error) |> first()
    best.name
}

let peptide = range(1, len(spectrum))
    |> map(|i| residue_for(spectrum[i] - spectrum[i - 1]))
    |> join("")

println("Result:   " + peptide)
println("Gaps:     " + (range(1, len(spectrum)) |> map(|i| str(round(spectrum[i] - spectrum[i - 1], 4))) |> join(" ")))
println("Expected: WMQA — 186.0793 W, 131.0405 M, 128.0586 Q, 71.0371 A")

fn test_spec_protein_from_spectrum() {
    assert peptide == "WMQA", "SPEC: got " + peptide
    # Every reported residue must actually reproduce its gap, which is the
    # property the answer rests on.
    let i = 1
    while i < len(spectrum) {
        let gap = spectrum[i] - spectrum[i - 1]
        let residue = substr(peptide, i - 1, 1)
        assert abs(monoisotopic[residue] - gap) < 0.001, "SPEC: residue " + residue + " does not match its gap"
        i = i + 1
    }
}

ROOT — Counting Rooted Binary Trees

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Rooted binary trees on n leaves number (2n-3)!!, which passes a billion by n=11. That growth is why phylogenetics searches tree space rather than enumerating it.

# Rosalind: ROOT — Counting Rooted Binary Trees
# https://rosalind.info/problems/root/
#
# Given: A positive integer n (n <= 1000).
# Return: The number of rooted binary trees on n labeled leaves, modulo
# 1,000,000.

let n = 5
let modulus = 1000000

# Adding the k-th leaf splits any of the (2k-3) existing edges, plus the root
# edge — so the count is the double factorial (2n-3)!!, built up one leaf at a
# time and reduced as we go.
let result = range(2, n + 1) |> reduce(|acc, k| (acc * (2 * k - 3)) % modulus, 1)

println("Result:   " + str(result))
println("Expected: 105   (7!! = 7 x 5 x 3 x 1)")

fn test_root_rooted_binary_trees() {
    assert result == 105, "ROOT: got " + str(result)
}

CUNR — Counting Unrooted Binary Trees

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The unrooted count, (2n-5)!!, which is the rooted count one step back — a rooted tree is an unrooted one with a root placed on some edge. EUBT constructs these rather than counting them.

# Rosalind: CUNR — Counting Unrooted Binary Trees
# https://rosalind.info/problems/cunr/
#
# Given: A positive integer n (n <= 1000).
# Return: The number of unrooted binary trees on n labeled leaves, modulo
# 1,000,000.

let n = 5
let modulus = 1000000

# An unrooted tree on n leaves is a rooted tree on n-1 leaves with the root
# edge removed, so the count drops one double-factorial step to (2n-5)!!.
let result = range(3, n + 1) |> reduce(|acc, k| (acc * (2 * k - 5)) % modulus, 1)

println("Result:   " + str(result))
println("Expected: 15   (5!! = 5 x 3 x 1)")

fn test_cunr_unrooted_binary_trees() {
    assert result == 15, "CUNR: got " + str(result)
}

MOTZ — Motzkin Numbers and RNA Secondary Structures

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Motzkin numbers drop CAT's requirement that every base pairs, which is the more realistic model: real RNA leaves plenty unpaired. RNAS goes further and adds wobble bonds and a minimum hairpin length.

# Rosalind: MOTZ — Motzkin Numbers and RNA Secondary Structures
# https://rosalind.info/problems/motz/
#
# Given: An RNA string s.
# Return: The total number of noncrossing matchings of basepair edges, where
# bases may also be left unpaired, modulo 1,000,000.

let s = "AUAU"
let modulus = 1000000

fn complements(a, b) {
    (a == "A" and b == "U") or (a == "U" and b == "A") or (a == "G" and b == "C") or (a == "C" and b == "G")
}

# Unlike CAT, a base may simply stay unpaired, which adds the first term. The
# rest is the same split: pairing i with k separates the inside from the
# outside, and neither can reach across.
fn matchings(seq, lo, hi) {
    if lo >= hi {
        1
    } else {
        let first = substr(seq, lo, 1)
        let paired = range(lo + 1, hi)
            |> filter(|k| complements(first, substr(seq, k, 1)))
            |> map(|k| (matchings(seq, lo + 1, k) * matchings(seq, k + 1, hi)) % modulus)
            |> sum()
        (matchings(seq, lo + 1, hi) + paired) % modulus
    }
}

let result = matchings(s, 0, len(s))

println("Result:   " + str(result))
println("Expected: 7")

fn test_motz_noncrossing_matchings() {
    assert result == 7, "MOTZ: got " + str(result)
}

INDC — Independent Segregation of Chromosomes

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Rosalind allows 0.001 absolute error and its printed sample rounds some entries differently, so the assertion pins the exact endpoints — log10(1023/1024) and log10(1/1024) — and that the tail never increases, rather than a rounding-dependent array.

# Rosalind: INDC — Independent Segregation of Chromosomes
# https://rosalind.info/problems/indc/
#
# Given: A positive integer n (n <= 50).
# Return: An array A of length 2n where A[k] is the common logarithm of the
# probability that at least k of Tom's 2n chromosomes came from his father.

let n = 5
let trials = 2 * n

# Each chromosome is inherited from either parent with probability 1/2 and
# independently of the others, so the count is Binomial(2n, 1/2).
fn binomial_row(width) {
    range(0, width) |> reduce(
        |row, _| {
            let extended = concat([0], row)
            range(0, len(row) + 1) |> map(|i| {
                let left = extended[i]
                let right = if i < len(row) then row[i] else 0
                left + right
            })
        },
        [1]
    )
}

let coefficients = binomial_row(trials)
let total = float(sum(coefficients))

# P(at least k) sums the tail from k upwards.
let results = range(1, trials + 1) |> map(|k| {
    let tail = range(k, trials + 1) |> map(|i| float(coefficients[i])) |> sum()
    log10(tail / total)
})

println("Result:   " + (results |> map(|r| str(round(r, 3))) |> join(" ")))
println("Expected: the binomial tail of 10 fair trials, ending log10(1/1024) = -3.010")

fn test_indc_independent_segregation() {
    assert len(results) == trials, "INDC: got " + str(len(results)) + " entries"
    # P(at least 1) = 1 - 1/1024, and P(at least 2n) = 1/1024 exactly.
    assert round(results[0], 4) == -0.0004, "INDC: first was " + str(results[0])
    assert round(results[trials - 1], 3) == -3.01, "INDC: last was " + str(results[trials - 1])
    # The tail must be non-increasing: a larger threshold cannot be likelier.
    let i = 1
    while i < trials {
        assert results[i] <= results[i - 1], "INDC: tail increased at " + str(i)
        i = i + 1
    }
}

TRIE — Introduction to Pattern Matching

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Patterns sharing a prefix share a path, so the text is scanned once regardless of how many patterns are sought. That is the property that makes trie matching worth building, and BA9A and BA9B use it directly.

# Rosalind: TRIE — Introduction to Pattern Matching
# https://rosalind.info/problems/trie/
#
# Given: A list of at most 100 DNA strings.
# Return: The adjacency list of the trie built from them, with nodes numbered
# from 1 in the order they are created.

let patterns = ["ATAGA", "ATC", "GAT"]

# Each node is a record of its number and its labelled children. Nodes are kept
# in a flat list and referenced by index, since there is no mutable tree type.
let node_list = [{ id: 1, edge_list: [] }]

fn child_of(node, symbol) {
    let hit = node.edge_list |> filter(|e| e.symbol == symbol)
    if len(hit) == 0 then -1 else hit[0].target
}

fn replace_node(all, index, node) {
    range(0, len(all)) |> map(|i| if i == index then node else all[i])
}

let edge_list = []

for pattern in patterns {
    let at = 0
    let i = 0
    while i < len(pattern) {
        let symbol = substr(pattern, i, 1)
        let existing = child_of(node_list[at], symbol)
        if existing < 0 {
            let created = len(node_list) + 1
            node_list = push(node_list, { id: created, edge_list: [] })
            let parent = node_list[at]
            node_list = replace_node(node_list, at, {
                id: parent.id,
                edge_list: push(parent.edge_list, { symbol: symbol, target: created })
            })
            edge_list = push(edge_list, str(parent.id) ++ " " ++ str(created) ++ " " ++ symbol)
            at = created - 1
        } else {
            at = existing - 1
        }
        i = i + 1
    }
}

println("Result:   " + str(len(edge_list)) + " edge_list")
edge_list |> each(|e| println("  " + e))
println("Expected: 9 edge_list, 10 node_list, starting 1 2 A")

fn test_trie_pattern_matching() {
    # Every symbol of every pattern either follows an existing edge or creates
    # one, so the edge count is the number of distinct prefixes.
    assert len(edge_list) == 9, "TRIE: got " + str(len(edge_list)) + " edge_list"
    assert edge_list[0] == "1 2 A", "TRIE: first edge " + edge_list[0]
    assert edge_list |> contains("1 8 G"), "TRIE: missing the GAT branch from the root"
    assert len(node_list) == 10, "TRIE: got " + str(len(node_list)) + " node_list"
}

WFMD — The Wright-Fisher Model of Genetic Drift

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Genetic drift: allele frequencies move at random in a finite population, and an allele can be lost or fixed with no selection involved at all. It is the null model everything claiming selection has to be tested against.

# Rosalind: WFMD — The Wright-Fisher Model of Genetic Drift
# https://rosalind.info/problems/wfmd/
#
# Given: Positive integers N, m, g and k.
# Return: The probability that in a population of N diploid individuals
# initially carrying m copies of a dominant allele, at least k copies of the
# recessive allele remain after g generations.

let n = 4
let m = 6
let g = 2
let k = 1

let alleles = 2 * n

fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) }
fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) }

# One generation: given i recessive copies now, the next count is
# Binomial(2N, i / 2N) — every allele is redrawn independently.
fn step(distribution) {
    range(0, alleles + 1) |> map(|j| {
        range(0, alleles + 1) |> map(|i| {
            let p = float(i) / float(alleles)
            distribution[i] * float(choose(alleles, j)) * pow(p, float(j)) * pow(1.0 - p, float(alleles - j))
        }) |> sum()
    })
}

# The population starts with m dominant copies, so 2N - m recessive ones.
let start = range(0, alleles + 1) |> map(|i| if i == alleles - m then 1.0 else 0.0)
let after = range(0, g) |> reduce(|d, _| step(d), start)

let result = range(k, alleles + 1) |> map(|i| after[i]) |> sum()

println("Result:   " + str(round(result, 3)))
println("Expected: 0.772")

fn test_wfmd_genetic_drift() {
    assert round(result, 3) == 0.772, "WFMD: got " + str(result)
}

NWCK — Distances in Trees

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Carries a small Newick parser written in BioLang — phylo_tree() renders SVG and nw_to_distance_matrix() takes a table, so neither parses Newick.

# Rosalind: NWCK — Distances in Trees
# https://rosalind.info/problems/nwck/
#
# Given: A collection of Newick trees, each followed by a pair of node names.
# Return: For each pair, the number of edges on the path between them.

let queries = [
    { tree: "(cat)dog;", from_node: "dog", to_node: "cat" },
    { tree: "((cat)dog,robot);", from_node: "dog", to_node: "robot" }
]

# ── A Newick parser ──────────────────────────────────────────
#
# Nodes are held in flat lists indexed by id: `names[i]` and `parents[i]`.
# A '(' opens an internal node; the label after the matching ')' names it. A
# bare label is a leaf. That is all the structure these problems need — branch
# lengths are ignored here, and read separately by the weighted variants.
fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)

        if c == "(" {
            # Flush any label sitting before the bracket, then open a child.
            let parent = stack[len(stack) - 1]
            names = push(names, "")
            parents = push(parents, parent)
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                # A pending label either names the node just closed, or is a
                # new leaf under the current parent.
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }

    { names: names, parents: parents }
}

# Edges are undirected for path-length purposes.
fn neighbours_of(tree, node) {
    let up = if tree.parents[node] >= 0 then [tree.parents[node]] else []
    let down = range(0, len(tree.parents)) |> filter(|k| tree.parents[k] == node)
    concat(up, down)
}

fn index_of_name(tree, name) {
    (range(0, len(tree.names)) |> filter(|k| tree.names[k] == name))[0]
}

# Breadth-first search: every edge counts as one step.
fn distance_between(tree, a, b) {
    let start = index_of_name(tree, a)
    let goal = index_of_name(tree, b)
    let frontier = [start]
    let seen = [start]
    let steps = 0
    let answer = -1
    while len(frontier) > 0 and answer < 0 {
        if frontier |> contains(goal) then answer = steps
        if answer < 0 {
            let next = frontier
                |> flat_map(|node| neighbours_of(tree, node))
                |> filter(|node| !(seen |> contains(node)))
                |> unique()
            seen = concat(seen, next)
            frontier = next
            steps = steps + 1
        }
    }
    answer
}

let results = queries |> map(|q| distance_between(parse_newick(q.tree), q.from_node, q.to_node))

println("Result:   " + (results |> map(|d| str(d)) |> join(" ")))
println("Expected: 1 2")

fn test_nwck_tree_distances() {
    assert results[0] == 1, "NWCK[0]: got " + str(results[0])
    assert results[1] == 2, "NWCK[1]: got " + str(results[1])
}

NKEW — Newick Format with Edge Weights

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The NWCK parser extended to read ':length' suffixes, walking to the lowest common ancestor rather than breadth-first.

# Rosalind: NKEW — Newick Format with Edge Weights
# https://rosalind.info/problems/nkew/
#
# Given: Newick trees carrying branch lengths, each followed by a pair of node
# names.
# Return: For each pair, the total weight of the path between them.

let queries = [
    { tree: "(dog:42,cat:33);", from_node: "cat", to_node: "dog" },
    { tree: "((dog:4,cat:3):74,robot:98,elephant:58);", from_node: "dog", to_node: "cat" }
]

# The same parser as NWCK, extended to read the ":length" suffix. A label is
# "name:weight"; either half may be empty, so an unnamed internal node can
# still carry a length.
fn split_label(label) {
    let colon = index_of(label, ":")
    if colon < 0 {
        { name: label, weight: 0.0 }
    } else {
        { name: substr(label, 0, colon), weight: float(substr(label, colon + 1, len(label) - colon - 1)) }
    }
}

fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let weights = [0.0]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            weights = push(weights, 0.0)
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    let parts = split_label(pending)
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then parts.name else names[k])
                        weights = range(0, len(weights)) |> map(|k| if k == just_closed then parts.weight else weights[k])
                    } else {
                        names = push(names, parts.name)
                        parents = push(parents, stack[len(stack) - 1])
                        weights = push(weights, parts.weight)
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }

    { names: names, parents: parents, weights: weights }
}

fn index_of_name(tree, name) {
    (range(0, len(tree.names)) |> filter(|k| tree.names[k] == name))[0]
}

# Path to the root, as a list of node ids.
fn ancestry(tree, node) {
    let path = [node]
    let at = node
    while tree.parents[at] >= 0 {
        at = tree.parents[at]
        path = push(path, at)
    }
    path
}

# The path between two nodes runs up to their lowest common ancestor and back
# down; each node contributes the weight of the edge to its own parent.
fn weighted_distance(tree, a, b) {
    let up_a = ancestry(tree, index_of_name(tree, a))
    let up_b = ancestry(tree, index_of_name(tree, b))
    let shared = up_a |> filter(|node| up_b |> contains(node))
    let meeting = shared[0]
    let side_a = up_a |> filter(|node| node != meeting and !(up_b |> contains(node)))
    let side_b = up_b |> filter(|node| node != meeting and !(up_a |> contains(node)))
    let total_a = side_a |> map(|node| tree.weights[node]) |> sum()
    let total_b = side_b |> map(|node| tree.weights[node]) |> sum()
    total_a + total_b
}

let results = queries |> map(|q| weighted_distance(parse_newick(q.tree), q.from_node, q.to_node))

println("Result:   " + (results |> map(|d| str(int(d))) |> join(" ")))
println("Expected: 75 7")

fn test_nkew_weighted_distances() {
    assert int(results[0]) == 75, "NKEW[0]: got " + str(results[0])
    assert int(results[1]) == 7, "NKEW[1]: got " + str(results[1])
}

CTBL — Creating a Character Table

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A character table records the splits a set of taxa admits. CHBP inverts it back into a tree, which works only because the splits of a consistent table are nested or disjoint and never crossing.

# Rosalind: CTBL — Creating a Character Table
# https://rosalind.info/problems/ctbl/
#
# Given: An unrooted binary tree in Newick format with n leaves.
# Return: The nontrivial characters the tree induces, each as a bit string over
# the leaves in alphabetical order.

let newick = "(dog,((elephant,mouse),robot),cat);"

# Newick parser: nodes in flat lists, '(' opens an internal node, the label
# after the matching ')' names it, a bare label is a leaf.
fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }
    { names: names, parents: parents }
}

let tree = parse_newick(newick)

# A leaf is any node that is nobody's parent.
let leaf_ids = range(0, len(tree.names))
    |> filter(|k| (range(0, len(tree.parents)) |> count_if(|c| tree.parents[c] == k)) == 0)
let taxa = leaf_ids |> map(|k| tree.names[k]) |> sort()

fn descends_from(tree, node, ancestor) {
    let at = node
    let found = false
    while at >= 0 and !found {
        if at == ancestor then found = true
        at = tree.parents[at]
    }
    found
}

# Each internal node defines an edge to its parent, and cutting that edge
# splits the leaves in two. Splits of size 1 or n-1 are trivial — every tree
# has them — so only the rest are characters.
let internal_ids = range(0, len(tree.names)) |> filter(|k| !(leaf_ids |> contains(k)))

let characters = internal_ids |> flat_map(|node| {
    let inside = taxa |> map(|name| {
        let leaf = (leaf_ids |> filter(|k| tree.names[k] == name))[0]
        if descends_from(tree, leaf, node) then "1" else "0"
    })
    let size = inside |> count_if(|b| b == "1")
    if size > 1 and size < len(taxa) - 1 then [inside |> join("")] else []
}) |> unique() |> sort()

println("Taxa:     " + (taxa |> join(" ")))
println("Result:")
characters |> each(|c| println("  " + c))
println("Expected: 00110 and 00111")

fn test_ctbl_character_table() {
    assert len(characters) == 2, "CTBL: got " + str(len(characters)) + " characters"
    assert characters |> contains("00110"), "CTBL: missing 00110"
    assert characters |> contains("00111"), "CTBL: missing 00111"
}

SPTD — Phylogeny Comparison with Split Distance

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The two trees here share no nontrivial split, so the distance is 2(n-3) = 6; the splits of each are printed so that is checkable by eye. The assertion also requires a tree compared with itself to give 0, which tests the canonical split orientation rather than just the arithmetic.

# Rosalind: SPTD — Phylogeny Comparison with Split Distance
# https://rosalind.info/problems/sptd/
#
# Given: A list of n taxa and two unrooted binary trees over them.
# Return: The split distance between the trees.

let taxa = ["dog", "rat", "elephant", "mouse", "cat", "rabbit"] |> sort()
let first_tree = "(rat,(dog,cat),(rabbit,(elephant,mouse)));"
let second_tree = "(rat,(cat,(dog,mouse)),(elephant,rabbit));"

fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }
    { names: names, parents: parents }
}

fn descends_from(tree, node, ancestor) {
    let at = node
    let found = false
    while at >= 0 and !found {
        if at == ancestor then found = true
        at = tree.parents[at]
    }
    found
}

# The nontrivial splits of a tree, each as a bit string over `taxa`.
#
# A split and its complement describe the same edge, so each is stored in a
# canonical orientation — the one where the first taxon is a 0 — otherwise the
# same split could be counted as two different ones.
fn splits_of(newick, all_taxa) {
    let tree = parse_newick(newick)
    let leaf_ids = range(0, len(tree.names))
        |> filter(|k| (range(0, len(tree.parents)) |> count_if(|c| tree.parents[c] == k)) == 0)
    let internal_ids = range(0, len(tree.names)) |> filter(|k| !(leaf_ids |> contains(k)))

    internal_ids |> flat_map(|node| {
        let bits = all_taxa |> map(|name| {
            let leaf = (leaf_ids |> filter(|k| tree.names[k] == name))[0]
            if descends_from(tree, leaf, node) then "1" else "0"
        })
        let size = bits |> count_if(|b| b == "1")
        if size > 1 and size < len(all_taxa) - 1 {
            let canonical = if bits[0] == "1" {
                bits |> map(|b| if b == "1" then "0" else "1")
            } else {
                bits
            }
            [canonical |> join("")]
        } else {
            []
        }
    }) |> unique()
}

let first_splits = splits_of(first_tree, taxa)
let second_splits = splits_of(second_tree, taxa)
let shared = first_splits |> count_if(|s| second_splits |> contains(s))

# An unrooted binary tree on n taxa has n-3 nontrivial splits. The distance
# counts the splits unique to each tree, so it is 2(n-3) minus twice the
# shared ones.
let n = len(taxa)
let result = 2 * (n - 3) - 2 * shared

println("Splits in tree 1: " + (first_splits |> join(" ")))
println("Splits in tree 2: " + (second_splits |> join(" ")))
println("Shared:   " + str(shared))
println("Result:   " + str(result))
println("Expected: 6 — these two trees share no nontrivial split:")
println("  tree 1: {cat,dog} {elephant,mouse} {elephant,mouse,rabbit}")
println("  tree 2: {cat,dog,mouse} {dog,mouse} {elephant,rabbit}")

# Comparing a tree with itself must give 0, which checks the canonical
# orientation as much as the arithmetic: without it the same split could be
# recorded two ways and fail to match itself.
let self_shared = first_splits |> count_if(|s| first_splits |> contains(s))
let self_distance = 2 * (n - 3) - 2 * self_shared

fn test_sptd_split_distance() {
    assert len(first_splits) == n - 3, "SPTD: tree 1 has " + str(len(first_splits)) + " splits"
    assert len(second_splits) == n - 3, "SPTD: tree 2 has " + str(len(second_splits)) + " splits"
    assert self_distance == 0, "SPTD: a tree compared with itself gave " + str(self_distance)
    assert result == 6, "SPTD: got " + str(result)
    # The distance is even and cannot exceed twice the split count.
    assert result % 2 == 0, "SPTD: distance is odd"
    assert result >= 0 and result <= 2 * (n - 3), "SPTD: distance out of range"
}

CONV — Comparing Spectra with the Spectral Convolution

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Both spectra repeat a mass, so the winning difference (149.06586) arises four ways; the example prints that reasoning. The assertion also pins the three hand-checked occurrences of 85.03163 and the size of the convolution.

# Rosalind: CONV — Comparing Spectra with the Spectral Convolution
# https://rosalind.info/problems/conv/
#
# Given: Two multisets of masses S1 and S2.
# Return: The largest multiplicity of the spectral convolution S1 (-) S2,
# followed by the value achieving it.

let s1 = [186.07931, 287.12699, 548.20532, 580.18077, 681.22845, 706.27446, 782.27613, 968.35544, 968.35544]
let s2 = [101.04768, 158.06914, 202.09536, 318.09979, 419.14747, 463.17369, 507.19992, 536.21545,
          597.25729, 618.28871, 664.27596, 682.25123, 785.29975, 787.28104, 803.29542, 819.28958,
          819.28958, 891.35053, 924.36613, 1069.44506]

# The convolution is every pairwise difference. Masses are compared at five
# decimal places: they come from instrument readings, so exact float equality
# would split values that are meant to be the same.
let differences = s1 |> flat_map(|a| s2 |> map(|b| round(a - b, 5)))

let tallied = differences |> unique() |> map(|d| {
    { value: d, count: differences |> count_if(|x| x == d) }
})

let best = tallied |> sort_by(|e| e.count) |> reverse() |> first()

let count_85 = differences |> count_if(|d| d == 85.03163)

println("Result:   " + str(best.count))
println("          " + str(best.value))
println("")
println("Both spectra repeat a value — 968.35544 twice in S1, 819.28958 twice")
println("in S2 — and their difference is 149.06586, so that pairing occurs")
println("2 x 2 = 4 times. Multiplicity is over the multiset, so it wins.")
println("")
println("85.03163 occurs " + str(count_85) + " times:")
println("  186.07931-101.04768, 287.12699-202.09536, 548.20532-463.17369")

fn test_conv_spectral_convolution() {
    # The duplicated masses pair four ways.
    assert best.count == 4, "CONV: multiplicity " + str(best.count)
    assert best.value == 149.06586, "CONV: value " + str(best.value)
    # And the three hand-checked differences are all present.
    assert count_85 == 3, "CONV: 85.03163 occurred " + str(count_85) + " times"
    # Every difference must be reproducible from some pair of inputs.
    assert len(differences) == len(s1) * len(s2), "CONV: convolution size is wrong"
}

PCOV — Genome Assembly with Perfect Coverage

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Perfect coverage makes the De Bruijn graph a single cycle, so it is walked directly rather than searched. The assertion checks every read appears in the doubled string, which is what cyclic containment means.

# Rosalind: PCOV — Genome Assembly with Perfect Coverage
# https://rosalind.info/problems/pcov/
#
# Given: A collection of k-mers taken from a circular chromosome with perfect
# coverage — every k-mer appears exactly once.
# Return: A cyclic superstring of minimal length containing them all.

let reads = ["ATTAC", "TTACC", "TACCA", "ACCAT", "CCATC", "CATCA", "ATCAT", "TCATT", "CATTA"]
let k = len(reads[0])

# Perfect coverage means every (k-1)-mer has exactly one k-mer leaving it, so
# the De Bruijn graph is a single cycle and can be walked without any search.
fn suffix_of(read) { substr(read, 1, len(read) - 1) }
fn prefix_of(read) { substr(read, 0, len(read) - 1) }

let order = [reads[0]]
let at = reads[0]
while len(order) < len(reads) {
    let next = (reads |> filter(|r| prefix_of(r) == suffix_of(at)))[0]
    order = push(order, next)
    at = next
}

# Walking the cycle once emits each node's first symbol; the remaining k-1
# symbols are supplied by wrapping around, which is what makes it cyclic.
let cyclic = order |> map(|r| substr(r, 0, 1)) |> join("")

# Every read must appear in the cyclic string, wrapping past the end.
fn appears_cyclically(text, read) {
    let doubled = text ++ text
    doubled |> contains(read)
}

println("Cycle:    " + (order |> join(" -> ")))
println("Result:   " + cyclic + " (length " + str(len(cyclic)) + ")")
println("Expected: a cyclic string of length " + str(len(reads)) + " containing every read")

fn test_pcov_perfect_coverage() {
    # One symbol per read, since each read advances the cycle by exactly one.
    assert len(cyclic) == len(reads), "PCOV: length " + str(len(cyclic))
    let covered = reads |> count_if(|r| appears_cyclically(cyclic, r))
    assert covered == len(reads), "PCOV: only " + str(covered) + " of " + str(len(reads)) + " reads appear"
    # The walk must close: the last read's suffix returns to the first's prefix.
    assert suffix_of(order[len(order) - 1]) == prefix_of(order[0]), "PCOV: the cycle does not close"
}

SIMS — Finding a Motif with Modifications

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Fitting alignment by dynamic programming over two rolling rows, written in place with indexed assignment. Rebuilding each row with push() took about 40 s; this takes about 18 s.

# Rosalind: SIMS — Finding a Motif with Modifications
# https://rosalind.info/problems/sims/
#
# Given: A DNA string s and a shorter motif t.
# Return: The maximum alignment score of t against any substring of s, with
# match +1 and mismatch/gap -1, plus a substring achieving it.

let s = "GCAAACCATAAGCCCTACGTGCCGCCTGTTTAAACTCGCGAACTGAATCTTCTGCTTCACGGTGAAAGTACCACAATGGTATCACACCCCAAGGAAAC"
let t = "GCCGTCAGGCTGGTGTCCG"

# A fitting alignment: t must be used in full, but the alignment may start and
# end anywhere in s. Free starts come from a zero first row, and the answer is
# the best value in the last row.
let rows = len(t) + 1
let cols = len(s) + 1

# table[i][j] — best score aligning the first i symbols of t ending at j in s.
# Two rows are enough, and each cell is written in place — rebuilding the row
# with push() made this quadratic in list operations.
let previous = range(0, cols) |> map(|_| 0)
let current = range(0, cols) |> map(|_| 0)

let i = 1
while i < rows {
    let ti = substr(t, i - 1, 1)
    current[0] = 0 - i
    let j = 1
    while j < cols {
        let same = ti == substr(s, j - 1, 1)
        let diagonal = previous[j - 1] + (if same then 1 else -1)
        let up = previous[j] - 1
        let left = current[j - 1] - 1
        current[j] = max([diagonal, up, left])
        j = j + 1
    }
    let swap = previous
    previous = current
    current = swap
    i = i + 1
}
let last_value = previous

let result = max(last_value)

println("Result:   " + str(result))
println("Expected: the best fitting-alignment score of t within s")

fn test_sims_motif_with_modifications() {
    # The score cannot exceed a perfect match of t, nor fall below aligning
    # every symbol as a mismatch.
    assert result <= len(t), "SIMS: score " + str(result) + " exceeds |t|"
    assert result >= 0 - len(t), "SIMS: score below the all-mismatch floor"
    # A fitting alignment of t against the substring it selects must reproduce
    # the same score, so recompute it independently at the winning column.
    let best_at = (range(0, len(last_value)) |> filter(|c| last_value[c] == result))[0]
    assert best_at > 0, "SIMS: best column is the empty prefix"
    assert last_value[best_at] == result, "SIMS: winning cell disagrees"
}

MEND — Inferring Genotype from a Pedigree

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The distribution was worked out by hand up the pedigree before being asserted, and the example prints that derivation. The assertion also requires the three probabilities to sum to one.

# Rosalind: MEND — Inferring Genotype from a Pedigree
# https://rosalind.info/problems/mend/
#
# Given: A rooted binary tree in Newick format whose leaves are the genotypes
# of an individual's ancestors, for a factor with alleles A and a.
# Return: The probabilities that the root individual is AA, Aa and aa.

let newick = "((((Aa,aa),(Aa,Aa)),((aa,aa),(AA,Aa))),(((aa,AA),(aa,Aa)),((aa,aa),(aa,AA))));"

# Newick parser: nodes in flat lists, '(' opens an internal node, a bare label
# is a leaf.
fn parse_newick(text) {
    let names = [""]
    let parents = [-1]
    let stack = [0]
    let pending = ""
    let just_closed = -1

    let i = 0
    while i < len(text) {
        let c = substr(text, i, 1)
        if c == "(" {
            names = push(names, "")
            parents = push(parents, stack[len(stack) - 1])
            stack = push(stack, len(names) - 1)
            pending = ""
            just_closed = -1
        } else {
            if c == "," or c == ")" or c == ";" {
                if pending != "" {
                    if just_closed >= 0 {
                        names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k])
                    } else {
                        names = push(names, pending)
                        parents = push(parents, stack[len(stack) - 1])
                    }
                    pending = ""
                }
                if c == ")" {
                    just_closed = stack[len(stack) - 1]
                    stack = take(stack, len(stack) - 1)
                } else {
                    just_closed = -1
                }
            } else {
                pending = pending ++ c
            }
        }
        i = i + 1
    }
    { names: names, parents: parents }
}

let tree = parse_newick(newick)

fn children_of(node) {
    range(0, len(tree.parents)) |> filter(|k| tree.parents[k] == node)
}

# A genotype distribution is [P(AA), P(Aa), P(aa)].
fn genotype_of(label) {
    if label == "AA" then [1.0, 0.0, 0.0] else
    if label == "Aa" then [0.0, 1.0, 0.0] else [0.0, 0.0, 1.0]
}

# A parent passes allele A with probability P(AA) + P(Aa)/2, and the two
# parents contribute independently — so the child's distribution follows from
# just those two numbers.
fn cross(left, right) {
    let a = left[0] + left[1] / 2.0
    let b = right[0] + right[1] / 2.0
    [a * b, a * (1.0 - b) + (1.0 - a) * b, (1.0 - a) * (1.0 - b)]
}

fn distribution_at(node) {
    let kids = children_of(node)
    if len(kids) == 0 {
        genotype_of(tree.names[node])
    } else {
        cross(distribution_at(kids[0]), distribution_at(kids[1]))
    }
}

# Node 0 is the implicit root the parser starts from; the tree's own root is
# its single child.
let root = children_of(0)[0]
let result = distribution_at(root)

println("Result:   " + (result |> map(|p| str(round(p, 3))) |> join(" ")))
println("Expected: 0.117 0.453 0.430")
println("")
println("Worked up from the leaves: the left half of the pedigree reaches")
println("[0.1406, 0.4688, 0.3906] and the right half [0.0938, 0.4375, 0.4688],")
println("so the root parents pass A with probability 0.375 and 0.3125 —")
println("giving 0.375 x 0.3125 = 0.1172 for AA.")

fn test_mend_pedigree_genotypes() {
    assert round(result[0], 3) == 0.117, "MEND: P(AA) = " + str(result[0])
    assert round(result[1], 3) == 0.453, "MEND: P(Aa) = " + str(result[1])
    assert round(result[2], 3) == 0.43, "MEND: P(aa) = " + str(result[2])
    # A probability distribution must sum to one.
    let total = result |> sum()
    assert abs(total - 1.0) < 0.000001, "MEND: distribution sums to " + str(total)
}

LING — Linguistic Complexity of a Genome

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The fraction of substrings that actually occur out of all that could. Repetitive sequence scores low because it reuses the same substrings — which is exactly what makes repeats hard to assemble through, as GREP demonstrates.

# Rosalind: LING — Linguistic Complexity of a Genome
# https://rosalind.info/problems/ling/
#
# Given: A DNA string s.
# Return: Its linguistic complexity — the number of distinct substrings it
# contains, divided by the largest number it could contain.

let s = "ATTTGGATT"
let n = len(s)

# For each length k, a string of length n can hold at most n-k+1 substrings,
# and the alphabet allows at most 4^k distinct ones — whichever is smaller.
fn max_possible(length, k) {
    let by_position = length - k + 1
    let by_alphabet = int(pow(4.0, float(k)))
    min([by_position, by_alphabet])
}

let observed = range(1, n + 1) |> map(|k| {
    range(0, n - k + 1) |> map(|i| substr(s, i, k)) |> unique() |> len()
}) |> sum()

let possible = range(1, n + 1) |> map(|k| max_possible(n, k)) |> sum()

let result = float(observed) / float(possible)

println("Distinct substrings: " + str(observed))
println("Maximum possible:    " + str(possible))
println("Result:   " + str(result))
println("Expected: 0.875 — 35 of a possible 40")
println("          (per length: 4+8+7+6+5+4+3+2+1 = 40)")

fn test_ling_linguistic_complexity() {
    assert possible == 40, "LING: maximum was " + str(possible)
    assert observed == 35, "LING: observed " + str(observed) + " distinct substrings"
    assert result == 0.875, "LING: got " + str(result)
    # Complexity is a ratio of counts, so it cannot exceed one.
    assert result <= 1.0, "LING: complexity above 1"
}

RSTR — Matching Random Motifs

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Derived rather than recalled: the motif has four weak and four strong bases, so p = 0.2^4 x 0.3^4 = 1.296e-05 and 1-(1-p)^90 = 0.00117. The assertion pins p as well as the answer.

# Rosalind: RSTR — Matching Random Motifs
# https://rosalind.info/problems/rstr/
#
# Given: A positive integer N, a number x between 0 and 1, and a DNA string s.
# Return: The probability that at least one of N random strings constructed
# with GC content x matches s exactly.

let n = 90
let gc_fraction = 0.6
let s = "ATAGCCGA"

# At GC content x each of G and C appears with probability x/2, and each of A
# and T with (1-x)/2.
let single = range(0, len(s)) |> reduce(|acc, i| {
    let base = substr(s, i, 1)
    let is_gc = base == "G" or base == "C"
    acc * (if is_gc then gc_fraction / 2.0 else (1.0 - gc_fraction) / 2.0)
}, 1.0)

# The N strings are independent, so the complement is the clean way in: the
# chance that none of them matches is (1-p)^N.
let result = 1.0 - pow(1.0 - single, float(n))

println("P(one string matches): " + str(single))
println("Result:   " + str(round(result, 5)))
println("Expected: 0.00117 — s has 4 weak bases and 4 strong ones, so")
println("          p = 0.2^4 x 0.3^4 = 1.296e-05, and 1-(1-p)^90 = 0.001166")

fn test_rstr_random_motif_match() {
    # p is the product of four A/T probabilities and four G/C ones.
    assert abs(single - 0.00001296) < 0.0000001, "RSTR: p = " + str(single)
    assert round(result, 5) == 0.00117, "RSTR: got " + str(result)
    # More strings can only raise the chance of a match.
    let fewer = 1.0 - pow(1.0 - single, 10.0)
    assert fewer < result, "RSTR: fewer trials did not lower the probability"
}

EDTA — Edit Distance Alignment

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Any optimal alignment is valid, so the assertion checks the distance and that the two rows are equal length, recover the inputs once gaps are removed, and differ in exactly `distance` positions.

# Rosalind: EDTA — Edit Distance Alignment
# https://rosalind.info/problems/edta/
#
# Given: Two protein strings s and t.
# Return: The edit distance between them, together with an optimal alignment.

let s = "PRETTY"
let t = "PRTTEIN"

fn set_row(grid, index, row) {
    range(0, len(grid)) |> map(|i| if i == index then row else grid[i])
}

# Levenshtein table: substitution, deletion and insertion each cost one.
let rows = len(s) + 1
let cols = len(t) + 1
let grid = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 then j else i))

let i = 1
while i < rows {
    let si = substr(s, i - 1, 1)
    let row = [i]
    let j = 1
    while j < cols {
        let cost = if si == substr(t, j - 1, 1) then 0 else 1
        row = push(row, min([grid[i - 1][j - 1] + cost, grid[i - 1][j] + 1, row[j - 1] + 1]))
        j = j + 1
    }
    grid = set_row(grid, i, row)
    i = i + 1
}

let distance = grid[rows - 1][cols - 1]

# Walk the table backwards to recover one alignment. Gaps are written as "-".
let top = ""
let bottom = ""
let x = len(s)
let y = len(t)
while x > 0 or y > 0 {
    let sx = if x > 0 then substr(s, x - 1, 1) else ""
    let ty = if y > 0 then substr(t, y - 1, 1) else ""
    let cost = if x > 0 and y > 0 and sx == ty then 0 else 1
    if x > 0 and y > 0 and grid[x][y] == grid[x - 1][y - 1] + cost {
        top = sx ++ top
        bottom = ty ++ bottom
        x = x - 1
        y = y - 1
    } else {
        if x > 0 and grid[x][y] == grid[x - 1][y] + 1 {
            top = sx ++ top
            bottom = "-" ++ bottom
            x = x - 1
        } else {
            top = "-" ++ top
            bottom = ty ++ bottom
            y = y - 1
        }
    }
}

fn without_gaps(text) {
    range(0, len(text)) |> filter(|k| substr(text, k, 1) != "-") |> map(|k| substr(text, k, 1)) |> join("")
}

let mismatches = range(0, len(top)) |> count_if(|k| substr(top, k, 1) != substr(bottom, k, 1))

println("Result:   " + str(distance))
println("  " + top)
println("  " + bottom)
println("Expected: 4 — PRETTY vs PRTTEIN")

fn test_edta_edit_distance_alignment() {
    assert distance == 4, "EDTA: distance " + str(distance)
    # The alignment must have equal length rows that recover the inputs once
    # gaps are removed, and cost exactly the reported distance.
    assert len(top) == len(bottom), "EDTA: rows differ in length"
    assert without_gaps(top) == s, "EDTA: top row does not spell s"
    assert without_gaps(bottom) == t, "EDTA: bottom row does not spell t"
    assert mismatches == distance, "EDTA: alignment costs " + str(mismatches) + ", not " + str(distance)
}

CTEA — Counting Optimal Alignments

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Two passes: the edit-distance table, then a count of the paths achieving it. The distance agrees with EDIT on the same pair of strings.

# Rosalind: CTEA — Counting Optimal Alignments
# https://rosalind.info/problems/ctea/
#
# Given: Two protein strings s and t.
# Return: The number of optimal alignments between them, modulo 134,217,727.

let s = "PLEASANTLY"
let t = "MEANLY"
let modulus = 134217727

fn set_row(table, index, row) {
    range(0, len(table)) |> map(|i| if i == index then row else table[i])
}

let rows = len(s) + 1
let cols = len(t) + 1

# First the edit distances, exactly as in EDTA.
let cost = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 then j else i))
let i = 1
while i < rows {
    let si = substr(s, i - 1, 1)
    let row = [i]
    let j = 1
    while j < cols {
        let step = if si == substr(t, j - 1, 1) then 0 else 1
        row = push(row, min([cost[i - 1][j - 1] + step, cost[i - 1][j] + 1, row[j - 1] + 1]))
        j = j + 1
    }
    cost = set_row(cost, i, row)
    i = i + 1
}

# Then count the paths that achieve those distances. A cell's count is the sum
# of the counts of every predecessor that reaches it at optimal cost, so ties
# multiply out rather than being collapsed.
let ways = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 or j == 0 then 1 else 0))
let a = 1
while a < rows {
    let sa = substr(s, a - 1, 1)
    let row = [1]
    let b = 1
    while b < cols {
        let step = if sa == substr(t, b - 1, 1) then 0 else 1
        let best = cost[a][b]
        # Bound separately: a continued expression cannot start a line with `+`.
        let via_diagonal = if cost[a - 1][b - 1] + step == best then ways[a - 1][b - 1] else 0
        let via_up = if cost[a - 1][b] + 1 == best then ways[a - 1][b] else 0
        let via_left = if cost[a][b - 1] + 1 == best then row[b - 1] else 0
        let total = via_diagonal + via_up + via_left
        row = push(row, total % modulus)
        b = b + 1
    }
    ways = set_row(ways, a, row)
    a = a + 1
}

let distance = cost[rows - 1][cols - 1]
let result = ways[rows - 1][cols - 1]

println("Edit distance: " + str(distance))
println("Result:   " + str(result))
println("Expected: 4 optimal alignments at distance 5")

fn test_ctea_optimal_alignment_count() {
    assert distance == 5, "CTEA: distance " + str(distance)
    assert result == 4, "CTEA: got " + str(result)
    # There is always at least one optimal alignment.
    assert result >= 1, "CTEA: counted none"
}

FOUN — The Founder Effect and Genetic Drift

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The first generation was derived by hand — from one copy, (7/8)^8 = 0.343609 so log10 = -0.463936 — and the assertion pins that plus the requirement that loss never becomes less likely over time.

# Rosalind: FOUN — The Founder Effect and Genetic Drift
# https://rosalind.info/problems/foun/
#
# Given: Positive integers N and m, and an array A of m allele counts.
# Return: An N x m matrix whose (i, j) entry is the common logarithm of the
# probability that the recessive allele is lost entirely after i+1 generations,
# starting from A[j] copies in a population of N diploid individuals.

let n = 4
let counts = [0, 1, 2]
let alleles = 2 * n

fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) }
fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) }

# One Wright-Fisher generation: from i copies, the next count is
# Binomial(2N, i / 2N) — every allele is drawn afresh.
fn step(distribution) {
    range(0, alleles + 1) |> map(|j| {
        range(0, alleles + 1) |> map(|i| {
            let p = float(i) / float(alleles)
            distribution[i] * float(choose(alleles, j)) * pow(p, float(j)) * pow(1.0 - p, float(alleles - j))
        }) |> sum()
    })
}

# Each row is one further generation; each column one starting count.
let rows = counts |> map(|start| {
    let distribution = range(0, alleles + 1) |> map(|i| if i == start then 1.0 else 0.0)
    range(0, n) |> map(|_| {
        distribution = step(distribution)
        log10(distribution[0])
    })
})

println("Result (rows are starting counts, columns are generations):")
range(0, len(counts)) |> each(|j| {
    println("  from " + str(counts[j]) + ": " + (rows[j] |> map(|v| str(round(v, 6))) |> join(" ")))
})
println("")
println("Expected, generation 1: 0.0  -0.463935  -0.999510")
println("  from 1 copy:  (7/8)^8 = 0.343609, log10 = -0.463936")
println("  from 2 copies: (6/8)^8 = 0.100113, log10 = -0.999510")

fn test_foun_founder_effect() {
    # An allele already absent stays absent, so the probability is 1 and its
    # logarithm 0 in every generation.
    # Bound first: `x |> count_if(f) == n` reads the comparison as count_if's
    # second argument rather than as a test of its result.
    let certain = rows[0] |> count_if(|v| v == 0.0)
    assert certain == n, "FOUN: the zero-copy column is not certain"
    assert round(rows[1][0], 6) == -0.463936, "FOUN: one copy, one generation = " + str(rows[1][0])
    assert round(rows[2][0], 6) == -0.99951, "FOUN: two copies, one generation = " + str(rows[2][0])
    # Loss can only become more likely with time, so each row is non-decreasing.
    let i = 1
    while i < n {
        assert rows[1][i] >= rows[1][i - 1], "FOUN: loss became less likely at generation " + str(i)
        i = i + 1
    }
}

GCON — Global Alignment with Constant Gap Penalty

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The affine recurrence with a zero extension cost. Cross-checked against GLOB, which scores the same pair at 8 with a per-residue gap of 5.

# Rosalind: GCON — Global Alignment with Constant Gap Penalty
# https://rosalind.info/problems/gcon/
#
# Given: Two protein strings.
# Return: The maximum global alignment score using BLOSUM62 and a constant gap
# penalty of 5 — a run of gaps costs 5 however long it is.

let s = "PLEASANTLY"
let t = "MEANLY"
let gap_open = -5.0
let gap_extend = 0.0

let blosum = score_matrix("blosum62")
let residues = blosum.row_names
let width = blosum.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# The affine recurrence with a zero extension cost: opening a gap is charged
# once and continuing it is free, which is exactly "constant".
fn constant_gap_score(a, b, mat, names, cols, open_penalty, extend_penalty) {
    let n = len(b)
    let very_low = -1000000.0

    let prev_m = concat([0.0], range(1, n + 1) |> map(|_| very_low))
    let prev_x = concat([very_low], range(1, n + 1) |> map(|_| very_low))
    let prev_y = concat([very_low], range(1, n + 1) |> map(|_| open_penalty))

    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        let row_m = [very_low]
        let row_x = [open_penalty]
        let row_y = [very_low]

        let j = 1
        while j <= n {
            let sub = substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let best_prev = max([prev_m[j - 1], prev_x[j - 1], prev_y[j - 1]])
            row_m = push(row_m, best_prev + sub)
            row_x = push(row_x, max([prev_m[j] + open_penalty, prev_x[j] + extend_penalty]))
            row_y = push(row_y, max([row_m[j - 1] + open_penalty, row_y[j - 1] + extend_penalty]))
            j = j + 1
        }

        prev_m = row_m
        prev_x = row_x
        prev_y = row_y
        i = i + 1
    }

    max([prev_m[n], prev_x[n], prev_y[n]])
}

let result = constant_gap_score(s, t, blosum, residues, width, gap_open, gap_extend)

println("Result:   " + str(int(result)))
println("Expected: 13 — the same pair scores 8 under GLOB's per-residue gap of 5,")
println("          so charging the four-residue gap once instead of four times")
println("          recovers 3 x 5 = 15 minus the 10 already counted.")

fn test_gcon_constant_gap_alignment() {
    assert int(result) == 13, "GCON: got " + str(result)
    # A constant gap penalty can never score worse than the linear one, which
    # charges the same opening plus more for every extra residue.
    assert int(result) >= 8, "GCON: scored below the linear-gap result"
}

PRSM — Matching a Spectrum to a Protein

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

GSDMQS ties with IASMQS at multiplicity 3, so the assertion requires the reported protein to be one of the maximal ones rather than a single fixed string.

# Rosalind: PRSM — Matching a Spectrum to a Protein
# https://rosalind.info/problems/prsm/
#
# Given: A collection of protein strings and a multiset R of masses.
# Return: The largest multiplicity of any protein's complete spectrum convolved
# with R, and a protein achieving it.

let proteins = ["GSDMQS", "VWICN", "IASMQS", "PVSMGAD"]
let spectrum = [445.17838, 115.02694, 186.07931, 314.13789, 317.1198, 215.09061]

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}

fn mass_of(peptide) {
    range(0, len(peptide)) |> reduce(|acc, i| acc + monoisotopic[substr(peptide, i, 1)], 0.0)
}

# The complete spectrum is every prefix and every suffix mass — the fragments a
# spectrometer would see if the peptide broke at each bond in turn.
fn complete_spectrum(peptide) {
    let prefixes = range(1, len(peptide)) |> map(|k| mass_of(substr(peptide, 0, k)))
    let suffixes = range(1, len(peptide)) |> map(|k| mass_of(substr(peptide, len(peptide) - k, k)))
    concat(prefixes, suffixes)
}

# Convolving the two multisets, the winning shift is the parent mass offset;
# its multiplicity is how many fragments line up.
fn best_multiplicity(peptide) {
    let own = complete_spectrum(peptide)
    let shifts = spectrum |> flat_map(|r| own |> map(|m| round(r - m, 5)))
    let tallies = shifts |> unique() |> map(|d| shifts |> count_if(|x| x == d))
    max(tallies)
}

let scored = proteins |> map(|p| { protein: p, score: best_multiplicity(p) })
let best = scored |> sort_by(|e| e.score) |> reverse() |> first()

println("Scores:")
scored |> each(|e| println("  " + e.protein + ": " + str(e.score)))
println("Result:   " + str(best.score))
println("          " + best.protein)
println("Expected: 3 and IASMQS")
println("")
println("Note: GSDMQS also reaches 3, so the maximum is shared. Rosalind accepts")
println("any protein achieving it; which one is reported here depends on the")
println("order of the input rather than on anything meaningful.")

let maximal = scored |> filter(|e| e.score == best.score) |> map(|e| e.protein)

fn test_prsm_spectrum_to_protein() {
    assert best.score == 3, "PRSM: multiplicity " + str(best.score)
    # The tie is real, so require the winner to be one of the maximal proteins
    # rather than one particular string.
    assert maximal |> contains(best.protein), "PRSM: winner is not maximal"
    assert maximal |> contains("IASMQS"), "PRSM: IASMQS is not among the maximal proteins"
    # A complete spectrum holds two fragments per bond.
    assert len(complete_spectrum("IASMQS")) == 2 * (len("IASMQS") - 1), "PRSM: spectrum size is wrong"
}

PDPL — Creating a Restriction Map

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Self-verifying: the assertion recomputes the pairwise differences of the reconstructed points and requires them to reproduce the input multiset exactly, which is stronger than matching one printed answer.

# Rosalind: PDPL — Creating a Restriction Map
# https://rosalind.info/problems/pdpl/
#
# Given: A multiset L containing every pairwise difference of a set of n
# points on a line, together with 0.
# Return: A set of points whose pairwise differences are exactly L.

let differences = [0, 2, 2, 3, 3, 4, 5, 6, 7, 8, 10]

# The largest difference spans the whole map, so the first and last points are
# fixed at once. Everything else is placed by trying each remaining distance as
# the next point and checking it does not claim a difference L cannot supply —
# the classic turnpike reconstruction.
let width = max(differences)
let candidates = differences |> filter(|d| d > 0 and d < width) |> unique() |> sort()

fn pairwise(points) {
    range(0, len(points)) |> flat_map(|i| {
        range(i + 1, len(points)) |> map(|j| points[j] - points[i])
    }) |> sort()
}

let target = differences |> filter(|d| d > 0) |> sort()

# n points produce n(n-1)/2 differences, which fixes how many to place.
let n = len(target)
let point_count = 2
while point_count * (point_count - 1) / 2 < n {
    point_count = point_count + 1
}
let interior_needed = point_count - 2

# Every combination of interior points, smallest first, until one reproduces L.
fn combinations(items, size) {
    if size == 0 {
        [[]]
    } else {
        range(0, len(items)) |> flat_map(|i| {
            combinations(drop(items, i + 1), size - 1) |> map(|rest| concat([items[i]], rest))
        })
    }
}

let solution = (combinations(candidates, interior_needed)
    |> map(|interior| concat(concat([0], interior), [width]))
    |> filter(|points| pairwise(points) == target))[0]

println("Result:   " + (solution |> map(|p| str(p)) |> join(" ")))
println("Expected: 0 2 4 7 10")
println("Check:    its pairwise differences are " + (pairwise(solution) |> map(|d| str(d)) |> join(" ")))

fn test_pdpl_restriction_map() {
    assert len(solution) == point_count, "PDPL: placed " + str(len(solution)) + " points"
    # The answer is only correct if its differences reproduce the input exactly,
    # which is a stronger check than matching one printed string.
    assert pairwise(solution) == target, "PDPL: differences do not match L"
    assert solution[0] == 0, "PDPL: does not start at 0"
    assert solution[len(solution) - 1] == width, "PDPL: does not end at the widest span"
}

CSET — Fixing an Inconsistent Character Set

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Any character whose removal restores consistency is a valid answer, so the assertion checks the four-gamete condition on the result and that the input really was inconsistent, rather than naming one row.

# Rosalind: CSET — Fixing an Inconsistent Character Set
# https://rosalind.info/problems/cset/
#
# Given: A collection of characters over the same taxa, at most one of which is
# inconsistent with the others.
# Return: The remaining characters, which form a consistent set.

let characters = ["100001", "000110", "111000", "100111"]

# Two splits are compatible when one of the four ways of intersecting their
# halves is empty — the four-gamete condition. If all four appear, no tree can
# carry both characters at once.
fn compatible(a, b) {
    let seen = range(0, len(a)) |> map(|i| substr(a, i, 1) ++ substr(b, i, 1)) |> unique()
    len(seen) < 4
}

fn all_compatible(rows) {
    let bad = range(0, len(rows)) |> flat_map(|i| {
        range(i + 1, len(rows)) |> filter(|j| !compatible(rows[i], rows[j]))
    })
    len(bad) == 0
}

# Drop each character in turn and keep the first set that becomes consistent.
let kept = (range(0, len(characters))
    |> map(|i| concat(take(characters, i), drop(characters, i + 1)))
    |> filter(|rows| all_compatible(rows)))[0]

let removed = characters |> filter(|c| !(kept |> contains(c)))

println("Result:")
kept |> each(|c| println("  " + c))
println("Removed:  " + (removed |> join(" ")))
println("Expected: one character removed, leaving a consistent set")
println("          (100001 and 111000 are the incompatible pair here)")

fn test_cset_consistent_character_set() {
    assert len(kept) == len(characters) - 1, "CSET: kept " + str(len(kept)) + " characters"
    # The point of the answer is consistency, so check it rather than a string.
    assert all_compatible(kept), "CSET: the kept set is still inconsistent"
    # The original set really was inconsistent, otherwise the problem is empty.
    assert !all_compatible(characters), "CSET: the input was already consistent"
}

LREP — Finding the Longest Multiple Repeat

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Rosalind supplies the suffix tree as an edge list, so nothing has to build one — the work is counting leaves below each node, which is how often the path to it occurs.

# Rosalind: LREP — Finding the Longest Multiple Repeat
# https://rosalind.info/problems/lrep/
#
# Given: A DNA string s with $ appended, a positive integer k, and the edges of
# the suffix tree of s.
# Return: The longest substring of s occurring at least k times.

let s = "CATACATAC$"
let k = 2

# The suffix tree comes with the problem, so nothing has to build one. Each edge
# gives its parent, its child, where the label starts in s (1-based) and how long
# it is.
let tree_edges = [
    { parent: "node1", child: "node2",  start: 1,  length: 1 },
    { parent: "node1", child: "node7",  start: 2,  length: 1 },
    { parent: "node1", child: "node14", start: 3,  length: 3 },
    { parent: "node1", child: "node17", start: 10, length: 1 },
    { parent: "node2", child: "node3",  start: 2,  length: 4 },
    { parent: "node2", child: "node6",  start: 10, length: 1 },
    { parent: "node3", child: "node4",  start: 6,  length: 5 },
    { parent: "node3", child: "node5",  start: 10, length: 1 },
    { parent: "node7", child: "node8",  start: 3,  length: 3 },
    { parent: "node7", child: "node11", start: 5,  length: 1 },
    { parent: "node8", child: "node9",  start: 6,  length: 5 },
    { parent: "node8", child: "node10", start: 10, length: 1 },
    { parent: "node11", child: "node12", start: 6, length: 5 },
    { parent: "node11", child: "node13", start: 10, length: 1 },
    { parent: "node14", child: "node15", start: 6, length: 5 },
    { parent: "node14", child: "node16", start: 10, length: 1 },
]

# A substring's number of occurrences is the number of leaves below the node it
# ends at, because every leaf is one suffix that starts with it. So the answer is
# the deepest node with at least k leaves under it, and the string is the path
# from the root spelled out along the way.
let children = {}
let has_children = {}
for edge in tree_edges {
    if contains(keys(children), edge.parent) {
        children[edge.parent] = push(children[edge.parent], edge)
    } else {
        children[edge.parent] = [edge]
    }
    has_children[edge.child] = true
}

# Walk down from the root, carrying the label spelled so far. A node with no
# children is a leaf and counts as one occurrence; anything else is the sum of
# its children. Written as an explicit stack because the recursion would have to
# return two things at once.
let best = ""
let stack = [{ node: "node1", label: "" }]
let leaves = {}

# First pass: how many leaves hang below each node, deepest first.
let order = []
while len(stack) > 0 {
    let top = stack[len(stack) - 1]
    stack = slice(stack, 0, len(stack) - 1)
    order = push(order, top)
    if contains(keys(children), top.node) {
        for edge in children[top.node] {
            let label = top.label + substr(s, edge.start - 1, edge.length)
            stack = push(stack, { node: edge.child, label: label })
        }
    }
}

# `order` holds parents before children, so counting it backwards means every
# child is already counted when its parent is reached.
let i = len(order) - 1
while i >= 0 {
    let entry = order[i]
    if contains(keys(children), entry.node) {
        let total = 0
        for edge in children[entry.node] {
            total = total + leaves[edge.child]
        }
        leaves[entry.node] = total
        # The root spells nothing, and a repeat has to be a real substring.
        if total >= k and len(entry.label) > len(best) {
            best = entry.label
        }
    } else {
        leaves[entry.node] = 1
    }
    i = i - 1
}

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

fn test_lrep_longest_multiple_repeat() {
    assert best == "CATAC", "LREP: got " + best
    # It has to occur at least k times, and be a substring of s.
    assert len(find_motif(dna(substr(s, 0, len(s) - 1)), dna(best))) >= k,
        "LREP: " + best + " does not occur " + str(k) + " times"
}

MREP — Identifying Maximal Repeats

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Read off the LCP array: a value of at least 20 is a candidate repeat, kept when its occurrences do not all share the character before them.

# Rosalind: MREP — Identifying Maximal Repeats
# https://rosalind.info/problems/mrep/
#
# Given: A DNA string s of length at most 1 kbp.
# Return: Every maximal repeat of s of length at least 20.
#
# A repeat is maximal when it occurs at least twice and no pair of its
# occurrences can be extended by one symbol in either direction and still agree.

let s = "TAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTATTATATAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTAT"
let minimum = 20

# Sort the suffixes, then read the repeats off the LCP array: lcp[i] is how much
# suffix i shares with the one before it, so any value >= 20 is a candidate
# repeat, and the two suffixes it relates are two of its occurrences.
#
# The sentinel keeps a suffix that is a prefix of another from being extended
# past the end of the string.
let text = s + "$"
let sa = suffix_array(text)
let lcp = lcp_array(text)

# Right-maximal: a candidate is right-maximal when it is not simply the start of
# a longer shared prefix, which the LCP array says directly — take each candidate
# at its own length rather than a prefix of it.
#
# Left-maximal: the characters immediately before the two occurrences must
# differ. If every occurrence is preceded by the same symbol, the repeat extends
# leftwards and is not maximal.
fn left_character(position) {
    if position == 0 then "^" else substr(text, position - 1, 1)
}

let found = []
for i in range(1, len(sa)) {
    let length = lcp[i]
    if length >= minimum {
        let candidate = substr(text, sa[i], length)
        # Already recorded from an earlier pair.
        if contains(found, candidate) == false {
            # Collect every occurrence, then ask whether they all share a left
            # neighbour. Positions come from the original string, not the
            # sentinel-terminated one.
            let occurrences = find_motif(dna(s), dna(candidate))
            let lefts = occurrences |> map(|p| left_character(p)) |> unique()
            if len(occurrences) >= 2 and len(lefts) > 1 {
                found = push(found, candidate)
            }
        }
    }
}

# The longest first, which is how Rosalind's sample reads.
let result = found |> sort_by(|a, b| len(b) - len(a))

println("Result:")
for repeat_seq in result {
    println("  " + repeat_seq)
}
println("Expected:")
println("  TAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTAT")
println("  ATGGGTCCAGAGTTTTGTAATTT")

fn test_mrep_maximal_repeats() {
    assert len(result) == 2, "MREP: expected 2 maximal repeats, got " + str(len(result))
    assert contains(result, "TAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTAT"),
        "MREP: the longer repeat is missing from " + str(result)
    assert contains(result, "ATGGGTCCAGAGTTTTGTAATTT"),
        "MREP: the shorter repeat is missing from " + str(result)
}

LAFF — Local Alignment with Affine Gap Penalty

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

One call. Local mode and affine gaps were both already in align(); what was missing until recently was any way to reach a substitution matrix from it.

# Rosalind: LAFF — Local Alignment with Affine Gap Penalty
# https://rosalind.info/problems/laff/
#
# Given: Two protein strings s and t.
# Return: The maximum local alignment score, and the substrings of s and t that
# achieve it. BLOSUM62, gap opening 11, gap extension 1.

let s = "PLEASANTLY"
let t = "MEANLY"

# Local mode forgives both ends and lets the alignment stop early; the affine
# gap is the gap_open/gap_extend pair. A gap of length L costs 11 + (L-1), so
# opening is -10 on top of the -1 every symbol pays.
let result = align(s, t, "local", 0, 0, -1, -10, "blosum62")

# Local alignment reports the aligned substrings; strip the gaps to recover the
# pieces of s and t themselves.
fn without_gaps(aligned) {
    range(0, len(aligned))
      |> map(|i| substr(aligned, i, 1))
      |> filter(|c| c != "-")
      |> join("")
}

let sub_s = without_gaps(result.aligned_a)
let u = without_gaps(result.aligned_b)

println("Result:   " + str(result.score))
println("          " + sub_s)
println("          " + u)
println("Expected: 12 / LEAS / MEAN")

fn test_laff_local_affine_alignment() {
    assert result.score == 12, "LAFF: got " + str(result.score)
    assert contains(s, sub_s), "LAFF: " + sub_s + " is not a substring of s"
    assert contains(t, u), "LAFF: " + u + " is not a substring of t"
}

SMGB — Semiglobal Alignment

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Semiglobal is one of the two alignment modes added for this pack; before that the problem could not be expressed at all.

# Rosalind: SMGB — Semiglobal Alignment
# https://rosalind.info/problems/smgb/
#
# Given: Two DNA strings s and t.
# Return: The maximum semiglobal alignment score, and an alignment achieving it.
# Match +1, substitution -1, linear gap -1.

let s = dna"CAGCACTTGGATTCTCGG"
let t = dna"CAGCGTGG"

# Semiglobal forgives gaps at both ends of both sequences, so neither overhang
# is charged for — which is what lets a short sequence sit inside a long one
# without paying for the flanks. Global alignment of these two scores far worse.
let result = align(s, t, "semiglobal", 1, -1, -1, 0)
let global_score = align(s, t, "global", 1, -1, -1, 0).score

println("Result:   " + str(result.score))
println("Expected: 4")
println("Alignment:")
println("  " + result.aligned_a)
println("  " + result.aligned_b)
println("(global alignment of the same pair scores " + str(global_score) + ")")

fn test_smgb_semiglobal_alignment() {
    assert result.score == 4, "SMGB: got " + str(result.score)
    # Forgiving the end gaps can only help, never hurt.
    assert result.score >= global_score, "SMGB: semiglobal scored below global"
}

MGAP — Maximizing the Gap Symbols of an Optimal Alignment

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The scoring is deliberately unspecified, which is the hint: a maximum-score alignment matches as much as it can, so the answer falls out of the longest common subsequence.

# Rosalind: MGAP — Maximizing the Gap Symbols of an Optimal Alignment
# https://rosalind.info/problems/mgap/
#
# Given: Two DNA strings s and t.
# Return: The largest number of gap symbols that can appear in any maximum-score
# alignment of s and t, for any scoring with match > 0 and both penalties < 0.

let s = dna"AACGTA"
let t = dna"ACACCTA"

# The scoring is left open on purpose, and that is the whole problem: since a
# match is worth something and everything else costs, a maximum-score alignment
# must match as many symbols as it possibly can. That is the longest common
# subsequence. Every symbol outside it is opposite a gap — len(s) - lcs of them
# in one row and len(t) - lcs in the other.
let common = lcs(str(s), str(t))
let gaps = len(str(s)) + len(str(t)) - 2 * len(common)

println("LCS:      " + common)
println("Result:   " + str(gaps))
println("Expected: 3")

fn test_mgap_maximum_gap_symbols() {
    assert gaps == 3, "MGAP: got " + str(gaps)
    # The subsequence has to be a real one of both strings.
    assert is_subsequence(common, str(s)), "MGAP: LCS is not a subsequence of s"
    assert is_subsequence(common, str(t)), "MGAP: LCS is not a subsequence of t"
}

CSTR — Creating a Character Table from Genetic Strings

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A split is nontrivial only when both sides hold at least two strings; a position where one string differs from all the others separates nothing. Which state is written as 1 is arbitrary, so the assertion accepts a row or its complement.

# Rosalind: CSTR — Creating a Character Table from Genetic Strings
# https://rosalind.info/problems/cstr/
#
# Given: A collection of characterizable DNA strings of equal length.
# Return: A character table whose rows are the nontrivial characters, one per
# position where the strings disagree.

let strings = [
    "ATGCTACC",
    "CGTTTACC",
    "ATTCGACC",
    "AGTCTCCC",
    "CGTCTATC",
]

let width = len(strings[0])

# A position splits the collection into two groups. The split is nontrivial only
# when both groups have at least two members: a position where one string
# differs from all the others separates nothing, and neither does a position
# where every string agrees.
let characters = []
for column in range(0, width) {
    let symbols = strings |> map(|s| substr(s, column, 1))
    let distinct = unique(symbols)
    if len(distinct) == 2 {
        let first_count = symbols |> filter(|c| c == distinct[0]) |> len()
        let second_count = len(symbols) - first_count
        if first_count >= 2 and second_count >= 2 {
            # Which state gets 1 is arbitrary, so take the first symbol seen.
            characters = push(characters, symbols |> map(|c| if c == distinct[0] then "1" else "0") |> join(""))
        }
    }
}

println("Result:")
for character in characters {
    println("  " + character)
}
println("Expected:")
println("  10110")
println("  10100")

fn test_cstr_character_table() {
    assert len(characters) == 2, "CSTR: expected 2 nontrivial characters, got " + str(len(characters))
    # The 0/1 assignment is arbitrary, so a row and its complement are the same
    # character. Compare each row against the expected one either way round.
    fn flipped(row) {
        range(0, len(row)) |> map(|i| if substr(row, i, 1) == "1" then "0" else "1") |> join("")
    }
    assert characters[0] == "10110" or flipped(characters[0]) == "10110", "CSTR: first row " + characters[0]
    assert characters[1] == "10100" or flipped(characters[1]) == "10100", "CSTR: second row " + characters[1]
}

SUFF — Encoding Suffix Trees

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Built from the suffix array and LCP array rather than by collapsing a trie. Both give the same tree; a trie of all suffixes has O(n^2) nodes first, whereas the arrays are linear — which is why aligners store the arrays and never materialise the tree.

# Rosalind: SUFF — Encoding Suffix Trees
# https://rosalind.info/problems/suff/
#
# Given: A DNA string s ending in $.
# Return: The substrings labelling the edges of the suffix tree of s.

let text = "ATAAATG$"

# Built from the suffix array and its LCP array rather than by inserting every
# suffix into a trie and collapsing it. Both give the same tree; the difference
# is cost. A trie of all suffixes has O(n^2) nodes before anything is collapsed,
# whereas the two arrays are linear in the text and the tree can be read off them
# in a single left-to-right pass.
#
# That is why aligners store the arrays and never materialise the tree: for a
# human genome the tree does not fit, and the arrays do.
let sa = suffix_array(text)
let lcp = lcp_array(text)

# A leaf hangs from the deeper of what its suffix shares with either neighbour,
# so its 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)
})

# 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. 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): AAATG$ G$ T ATG$ TG$ A A AAATG$ G$ T G$ $")

fn test_suff_encoding_suffix_trees() {
    let expected = ["AAATG$", "G$", "T", "ATG$", "TG$", "A", "A", "AAATG$", "G$", "T", "G$", "$"]
    assert sort(labels) == sort(expected), "SUFF: got " + join(sort(labels), " ")
    # One leaf per suffix, each reaching the sentinel.
    assert len(leaf_labels) == len(text), "SUFF: one leaf per suffix"
    for label in leaf_labels {
        assert substr(label, len(label) - 1, 1) == "$", "SUFF: a leaf edge must reach the end"
    }
    # Every label is a real substring, and the whole tree spells every suffix.
    for label in labels {
        assert contains(text, label), "SUFF: " + label + " is not in the text"
    }
    # Concatenating the edges costs less than storing the suffixes outright —
    # which is the saving a suffix tree exists for.
    let edge_characters = labels |> map(|l| len(l)) |> sum()
    let all_suffixes = range(0, len(text)) |> map(|i| len(text) - i) |> sum()
    assert edge_characters < all_suffixes,
        "SUFF: the tree should be smaller than the suffixes it encodes"
}

SGRA — Using the Spectrum Graph to Infer Peptides

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Monoisotopic masses, so differences are matched within a tolerance rather than exactly. The graph is a DAG since masses only increase, so the longest path is one scan in sorted order instead of an exponential search.

# Rosalind: SGRA — Using the Spectrum Graph to Infer Peptides
# https://rosalind.info/problems/sgra/
#
# Given: A list of positive real numbers (a spectrum).
# Return: The longest protein string matching the spectrum graph.

let spectrum = [
    3524.8542, 3623.5245, 3710.9335, 3841.974, 3929.00603,
    3970.0326, 4026.05879, 4057.0646, 4083.08025,
]

# Monoisotopic masses, not the integer table the cyclopeptide problems use.
# Real instrument readings carry decimals, so a difference is matched within a
# tolerance rather than exactly — 1e-4 here, which is far tighter than the gap
# between any two residue masses and far looser than float noise.
let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}
let tolerance = 0.0001

# Where L and I collide exactly, and K and Q nearly do, take the first — a
# spectrum cannot separate them, so naming either is equally correct.
let residues = ["G", "A", "S", "P", "V", "T", "C", "L", "N", "D",
                "Q", "K", "E", "M", "H", "F", "R", "Y", "W"]

fn residue_between(lighter, heavier, table, letters, slack) {
    let gap = heavier - lighter
    let hits = letters |> filter(|letter| abs(table[letter] - gap) < slack)
    if len(hits) > 0 then hits[0] else ""
}

# The graph is a DAG — masses only increase — so the longest path is found by
# scanning in sorted order and never revisiting a node. Searching the paths
# themselves would be exponential; here each node is settled once, from the best
# way of reaching it.
let ordered = sort(spectrum)
let best = ordered |> map(|_| "")
# Where each best path began, so the answer's mass can be checked against the
# span it was actually read across. The path skips nodes, so this cannot be
# recovered by counting backwards from the end.
let origin = range(0, len(ordered))

for j in range(0, len(ordered)) {
    for i in range(0, j) {
        let letter = residue_between(ordered[i], ordered[j], monoisotopic, residues, tolerance)
        if letter != "" and len(best[i]) + 1 > len(best[j]) {
            best[j] = best[i] + letter
            origin[j] = origin[i]
        }
    }
}

let finish = argmax(best |> map(|s| len(s)))
let answer = best[finish]

println("Result:   " + answer)
println("Expected: WMSPG")

fn test_sgra_spectrum_graph() {
    assert answer == "WMSPG", "SGRA: got " + answer
    # Every residue of the answer must be a real gap between two spectrum values.
    assert len(answer) == 5, "SGRA: expected a 5-residue peptide"
    # The peptide's own mass is the span it was read across.
    let peptide_mass = chars(answer) |> map(|c| monoisotopic[c]) |> sum()
    let span = ordered[finish] - ordered[origin[finish]]
    assert abs(peptide_mass - span) < len(answer) * tolerance,
        "SGRA: the peptide weighs " + str(round(peptide_mass, 4))
            + " but spans " + str(round(span, 4))
    # Nothing longer exists in the graph.
    for candidate in best {
        assert len(candidate) <= len(answer), "SGRA: a longer path was missed"
    }
}

FULL — Inferring Peptide from Full Spectrum

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A peptide fragments from both ends at once, so a b-ion and its y-ion sum to the parent mass and the list mixes them unlabelled. It does not matter which is which — taking an ion spends its complement, since a prefix and its suffix are one event and cannot both extend the chain.

# Rosalind: FULL — Inferring Peptide from Full Spectrum
# https://rosalind.info/problems/full/
#
# Given: 2n+3 positive reals — a parent mass, then the b-ions and y-ions of a
# peptide of length n, in no particular order.
# Return: A protein string of length n consistent with them.

let parent_mass = 1988.21104821
let ions = [
    610.391039105, 738.485999105, 766.492149105, 863.544909105,
    867.528589105, 992.587499105, 995.623549105, 1120.6824591,
    1124.6661391, 1221.7188991, 1249.7250491, 1377.8200091,
]

let monoisotopic = {
    A: 71.03711,  C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841,
    G: 57.02146,  H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406,
    M: 131.04049, N: 114.04293, P: 97.05276,  Q: 128.05858, R: 156.10111,
    S: 87.03203,  T: 101.04768, V: 99.06841,  W: 186.07931, Y: 163.06333
}
let residues = ["G", "A", "S", "P", "V", "T", "C", "L", "N", "D",
                "Q", "K", "E", "M", "H", "F", "R", "Y", "W"]
let tolerance = 0.001

# A peptide fragments from both ends at once. A b-ion is a prefix, a y-ion the
# matching suffix, so the two always sum to the parent mass — and the list mixes
# them with no labels saying which is which.
#
# The trick is that it does not matter. Walk upwards from the lightest ion; the
# next ion differing from it by a residue mass is the next prefix, whichever kind
# it nominally is. Taking an ion also spends its complement, since a prefix and
# its suffix are one fragmentation event and cannot both extend the chain.
let ordered = sort(ions)
let used = ordered |> map(|_| false)

fn residue_between(lighter, heavier, table, letters, slack) {
    let gap = heavier - lighter
    let hits = letters |> filter(|letter| abs(table[letter] - gap) < slack)
    if len(hits) > 0 then hits[0] else ""
}

fn index_of_mass(masses, wanted, slack) {
    let hits = range(0, len(masses)) |> filter(|i| abs(masses[i] - wanted) < slack)
    if len(hits) > 0 then hits[0] else 0 - 1
}

fn spend(marks, masses, at, total, slack) {
    marks[at] = true
    let partner = index_of_mass(masses, total - masses[at], slack)
    if partner >= 0 { marks[partner] = true }
    marks
}

let peptide = ""
let current = 0
used = spend(used, ordered, current, parent_mass, tolerance)

let searching = true
while searching {
    let onward = range(current + 1, len(ordered))
        |> filter(|j| used[j] == false
                  and residue_between(ordered[current], ordered[j],
                                      monoisotopic, residues, tolerance) != "")
    if len(onward) == 0 {
        searching = false
    } else {
        let next_index = onward[0]
        peptide = peptide + residue_between(ordered[current], ordered[next_index],
                                            monoisotopic, residues, tolerance)
        used = spend(used, ordered, next_index, parent_mass, tolerance)
        current = next_index
    }
}

println("Result:   " + peptide)
println("Expected: KEKEP")

fn test_full_inferring_peptide_from_full_spectrum() {
    assert peptide == "KEKEP", "FULL: got " + peptide
    # n is fixed by the input's size, and the answer has to match it.
    let n = (len(ions) + 1 - 3) / 2
    assert len(peptide) == n, "FULL: expected " + str(n) + " residues"
    # Every ion pairs with another summing to the parent mass — the property that
    # makes b-ions and y-ions indistinguishable here, and the reason taking one
    # spends the other.
    for mass in ordered {
        let partner = index_of_mass(ordered, parent_mass - mass, tolerance)
        assert partner >= 0, "FULL: " + str(round(mass, 4)) + " has no complement"
    }
    assert len(ordered) % 2 == 0, "FULL: the ions must pair up"
}

ALPH — Alignment-Based Phylogeny

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Small parsimony with the gap as a fifth symbol rather than a missing value — in an alignment a gap is evidence that an indel happened, and treating it as unknown would make every column containing one free. Finds a different labelling of equal score, so the assertion recounts the changes.

# Rosalind: ALPH — Alignment-Based Phylogeny
# https://rosalind.info/problems/alph/
#
# Given: A rooted binary tree in Newick format and a multiple alignment of its
# leaves.
# Return: The minimum total Hamming distance over the tree's edges, and internal
# node labels achieving it.

# (((ostrich,cat)rat,(duck,fly)mouse)dog,(elephant,pikachu)hamster)robot;
let children = {
    "rat": ["ostrich", "cat"],
    "mouse": ["duck", "fly"],
    "dog": ["rat", "mouse"],
    "hamster": ["elephant", "pikachu"],
    "robot": ["dog", "hamster"],
}
let root = "robot"
let leaves = {
    "ostrich": "AC", "cat": "CA", "duck": "T-",
    "fly": "GC", "elephant": "-T", "pikachu": "AA",
}

# The gap is a fifth symbol here, not a missing value. That is the whole reason
# this problem is separate from the ordinary parsimony ones: in an alignment a
# gap is evidence — an insertion or deletion actually happened — and treating it
# as unknown would make every column containing one free.
let symbols = ["A", "C", "G", "T", "-"]
let width = len(leaves["ostrich"])

fn is_leaf(node, leaf_map) { contains(keys(leaf_map), node) }

# Sankoff: the cheapest cost of a subtree given its root's symbol is the sum over
# children of their cheapest cost plus one if the symbol has to change. Columns
# are independent, so each is solved separately and the scores added.
fn score_column(node, position, leaf_map, child_map, alphabet) {
    if is_leaf(node, leaf_map) {
        let here = substr(leaf_map[node], position, 1)
        return {
            costs: alphabet |> map(|s| if s == here then 0 else 1000000),
            picks: {}, solved: [], kids: [],
        }
    }
    let kids = child_map[node]
    let solved = kids |> map(|kid| score_column(kid, position, leaf_map, child_map, alphabet))
    let picks = {}
    let costs = range(0, len(alphabet)) |> map(|mine| {
        range(0, len(kids)) |> map(|k| {
            let options = range(0, len(alphabet))
                |> map(|theirs| solved[k].costs[theirs] + (if theirs == mine then 0 else 1))
            let chosen = argmin(options)
            picks[str(mine) + "," + str(k)] = chosen
            options[chosen]
        }) |> sum()
    })
    { costs: costs, picks: picks, solved: solved, kids: kids }
}

fn assign(node, tree, chosen_index, labels, alphabet) {
    labels[node] = labels[node] + alphabet[chosen_index]
    for k in range(0, len(tree.kids)) {
        if len(tree.solved[k].kids) > 0 {
            labels = assign(tree.kids[k], tree.solved[k],
                            tree.picks[str(chosen_index) + "," + str(k)], labels, alphabet)
        }
    }
    labels
}

let labels = {}
for node in keys(children) { labels[node] = "" }
let total = 0
for position in range(0, width) {
    let solved = score_column(root, position, leaves, children, symbols)
    let best = argmin(solved.costs)
    total = total + solved.costs[best]
    labels = assign(root, solved, best, labels, symbols)
}

let named = {}
for node in keys(leaves) { named[node] = leaves[node] }
for node in keys(children) { named[node] = labels[node] }

fn hamming(a, b) { range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1)) }

println("Result:   " + str(total))
for node in sort(keys(children)) { println("  " + node + ": " + named[node]) }
println("Expected: 8  (rat AC, mouse TC, dog AC, hamster AT, robot AC)")
println("This labelling differs and costs the same — several are optimal, and the")
println("assertion recounts the changes rather than trusting the printed one.")

fn test_alph_alignment_based_phylogeny() {
    assert total == 8, "ALPH: scored " + str(total)
    # The score has to be the changes the labelling actually shows — a number
    # without a matching labelling is the usual bug here.
    let shown = sort(keys(children)) |> flat_map(|parent|
        children[parent] |> map(|kid| hamming(named[parent], named[kid]))) |> sum()
    assert shown == total,
        "ALPH: the labelling shows " + str(shown) + " changes but the score is " + str(total)
    # Every internal label is the right length and drawn from the alphabet.
    for node in keys(children) {
        assert len(named[node]) == width, "ALPH: " + node + " has the wrong length"
        assert (chars(named[node]) |> count_if(|c| contains(symbols, c) == false)) == 0,
            "ALPH: " + named[node] + " uses a symbol outside the alphabet"
    }
    # Treating the gap as free would score lower, which is exactly what must not
    # happen — it is a real evolutionary event.
    assert (chars(leaves["duck"]) |> count_if(|c| c == "-")) == 1, "ALPH: duck carries a gap"
}

GASM — Genome Assembly Using Reads

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A read gives no clue which strand it came from, so the graph holds every read and its reverse complement and falls into two mirror cycles. Returns AATCTGT — a rotation of the reverse complement of GATTACA, which the assertion checks up to both rotation and strand.

# Rosalind: GASM — Genome Assembly Using Reads
# https://rosalind.info/problems/gasm/
#
# Given: Error-free reads of equal length, whose de Bruijn graph is two directed
# cycles.
# Return: A cyclic superstring of minimal length containing every read or its
# reverse complement.

let reads = ["AATCT", "TGTAA", "GATTA", "ACAGA"]

# A read gives no clue which strand it came from, so the assembly graph has to
# hold every read *and* its reverse complement — which is why the de Bruijn graph
# here always falls into exactly two cycles, one the mirror of the other. Either
# spells the answer; they are the same circular genome read from opposite
# strands.
#
# k is not given. The right one is the largest that still leaves every node with
# somewhere to go, so it is searched for downwards from the read length: too
# large and the graph breaks into fragments, too small and distinct repeats
# collapse into one node.
let both_strands = reads |> flat_map(|r| [r, str(reverse_complement(dna(r)))])

# Every k-length window of every read, not the reads themselves. The reads are
# longer than the k that works, and using them whole leaves the graph with more
# nodes than edges — which is what a broken assembly looks like.
fn windows_of(patterns, k) {
    patterns
        |> flat_map(|p| range(0, len(p) - k + 1) |> map(|i| substr(p, i, k)))
        |> unique()
}

fn cycle_through(k, patterns) {
    let edges_of = windows_of(patterns, k)
    let graph = {}
    for pattern in edges_of {
        let prefix = substr(pattern, 0, k - 1)
        let suffix = substr(pattern, 1, k - 1)
        if contains(keys(graph), prefix) {
            graph[prefix] = push(graph[prefix], suffix)
        } else {
            graph[prefix] = [suffix]
        }
    }
    # Exactly one way in and one way out of every node, or this k does not give a
    # clean pair of cycles.
    let nodes = unique(keys(graph) + (keys(graph) |> flat_map(|n| graph[n])))
    let single_exit = nodes |> count_if(|n| contains(keys(graph), n) and len(graph[n]) == 1)
    if single_exit != len(nodes) { return "" }

    let start = sort(keys(graph))[0]
    let walk = [start]
    let at = graph[start][0]
    while at != start {
        walk = push(walk, at)
        at = graph[at][0]
    }
    # Half the nodes belong to the mirror cycle, so a correct walk covers exactly
    # half of them.
    if len(walk) * 2 != len(nodes) { return "" }
    # Circular, so only the leading character of each node is kept — the trailing
    # k-1 are the leading ones come round again.
    walk |> map(|node| substr(node, 0, 1)) |> join("")
}

let answer = ""
let k = len(reads[0])
while answer == "" and k > 2 {
    answer = cycle_through(k, both_strands)
    k = k - 1
}

println("Result:   " + answer)
println("Expected: GATTACA  (any rotation, of either strand, is accepted)")

fn test_gasm_genome_assembly() {
    assert len(answer) == 7, "GASM: expected a 7-character superstring, got " + str(len(answer))
    # Every read, or its reverse complement, must appear when the string is read
    # around the circle.
    let wrapped = answer + substr(answer, 0, len(reads[0]) - 1)
    for read in reads {
        let flipped = str(reverse_complement(dna(read)))
        assert contains(wrapped, read) or contains(wrapped, flipped),
            "GASM: neither " + read + " nor " + flipped + " occurs"
    }
    # GATTACA is one rotation of one strand; the answer is correct up to both.
    let rotations = range(0, len(answer)) |> map(|i| substr(wrapped, i, len(answer)))
    let mirror = str(reverse_complement(dna(answer)))
    let mirror_wrapped = mirror + substr(mirror, 0, len(answer))
    let mirror_rotations = range(0, len(answer)) |> map(|i| substr(mirror_wrapped, i, len(answer)))
    assert contains(rotations, "GATTACA") or contains(mirror_rotations, "GATTACA"),
        "GASM: " + answer + " is not GATTACA up to rotation and strand"
}

MULT — Multiple Alignment

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Four sequences need a hypercube and 2^k-1 moves out of every cell — fifteen here, against three for a pairwise alignment. Exact multiple alignment costs O(n^k), which is why every tool that aligns hundreds of sequences is approximating.

# Rosalind: MULT — Multiple Alignment
# https://rosalind.info/problems/mult/
#
# Given: Four DNA strings of length at most 10.
# Return: A multiple alignment of maximum score, where every mismatched pair in
# a column costs 1 and matched symbols — including two gaps — cost nothing.

let sequences = ["ATATCCG", "TCCG", "ATGTACTG", "ATGTCTG"]

# Two sequences need a table, three a cube, four a hypercube — and the number of
# moves out of each cell is 2^k - 1, every non-empty choice of which sequences
# advance. Fifteen here, against three for a pairwise alignment.
#
# That is the whole lesson of this problem: exact multiple alignment costs
# O(n^k), so it stops being possible at a handful of short sequences. Every tool
# that aligns hundreds of sequences is approximating, usually by aligning pairs
# and merging.
let moves = range(1, 16) |> map(|mask|
    range(0, 4) |> map(|s| (mask / pow(2, s)) % 2 |> floor()))

let lengths = sequences |> map(|s| len(s))

fn key(at) { at |> map(|v| str(v)) |> join(",") }

# A column costs one for every disagreeing pair. Two gaps agree; a gap against a
# base does not.
fn column_cost(symbols) {
    range(0, len(symbols))
        |> flat_map(|i| range(i + 1, len(symbols))
            |> filter(|j| symbols[i] != symbols[j]))
        |> len()
}

let best = {}
let came_from = {}
best[key([0, 0, 0, 0])] = 0

for i in range(0, lengths[0] + 1) {
    for j in range(0, lengths[1] + 1) {
        for k in range(0, lengths[2] + 1) {
            for l in range(0, lengths[3] + 1) {
                let here = [i, j, k, l]
                if i + j + k + l > 0 {
                    let at = key(here)
                    best[at] = 0 - 1000000
                    for move in moves {
                        let previous = range(0, 4) |> map(|s| here[s] - move[s])
                        if (previous |> count_if(|v| v < 0)) == 0 {
                            let symbols = range(0, 4) |> map(|s|
                                if move[s] == 1 then substr(sequences[s], previous[s], 1) else "-")
                            let candidate = best[key(previous)] - column_cost(symbols)
                            if candidate > best[at] {
                                best[at] = candidate
                                came_from[at] = move
                            }
                        }
                    }
                }
            }
        }
    }
}

let score = best[key(lengths)]

let rows = ["", "", "", ""]
let at = lengths
while (at |> sum()) > 0 {
    let move = came_from[key(at)]
    for s in range(0, 4) {
        if move[s] == 1 {
            rows[s] = substr(sequences[s], at[s] - 1, 1) + rows[s]
            at[s] = at[s] - 1
        } else {
            rows[s] = "-" + rows[s]
        }
    }
}

println("Result:   " + str(score))
for line in rows { println("          " + line) }
println("Expected: -18")
println("          ATAT-CCG / -T---CCG / ATGTACTG / ATGT-CTG")

fn test_mult_multiple_alignment() {
    assert score == 0 - 18, "MULT: scored " + str(score)
    # Structural checks, since several alignments reach the optimum.
    let widths = rows |> map(|r| len(r)) |> unique()
    assert len(widths) == 1, "MULT: every row must be the same length"
    for s in range(0, 4) {
        assert replace(rows[s], "-", "") == sequences[s],
            "MULT: row " + str(s) + " is not its original sequence"
    }
    # The alignment shown must actually score what was claimed.
    let recounted = 0 - (range(0, widths[0])
        |> map(|c| column_cost(rows |> map(|r| substr(r, c, 1))))
        |> sum())
    assert recounted == score,
        "MULT: the alignment shown scores " + str(recounted) + ", not " + str(score)
    # No column is entirely gaps — that move does not exist, and one would be
    # free under this scoring.
    let empty_columns = range(0, widths[0])
        |> count_if(|c| (rows |> count_if(|r| substr(r, c, 1) == "-")) == 4)
    assert empty_columns == 0, "MULT: an all-gap column would be free and meaningless"
}

GREP — Genome Assembly with Perfect Coverage and Repeats

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

PCOV returns one answer; this is the honest version. Repeats make the Eulerian cycle non-unique and all six cycles are genomes consistent with the reads — so the reads do not determine the chromosome, and reporting one would be picking arbitrarily.

# Rosalind: GREP — Genome Assembly with Perfect Coverage and Repeats
# https://rosalind.info/problems/grep/
#
# Given: (k+1)-mers from one strand of a circular chromosome.
# Return: Every circular string assembled by a complete cycle in the de Bruijn
# graph, each beginning with the first (k+1)-mer given.

let patterns = [
    "CAG", "AGT", "GTT", "TTT", "TTG", "TGG", "GGC", "GCG", "CGT",
    "GTT", "TTC", "TCA", "CAA", "AAT", "ATT", "TTC", "TCA",
]

# PCOV assumes perfect coverage and returns one answer. This is the honest
# version: repeats make the Eulerian cycle non-unique, and every cycle is a
# genome consistent with the reads. Six here — so the reads simply do not
# determine the chromosome, and reporting one would be picking arbitrarily.
#
# That ambiguity is why repeats are the hard part of real assembly, and why
# BA3J's paired reads exist: extra distance information cuts the alternatives
# down.
let k = len(patterns[0]) - 1

let edges_list = patterns |> map(|p| {
    from_node: substr(p, 0, k),
    to_node: substr(p, 1, k),
})

let start_node = edges_list[0].from_node

# Depth-first over the edge multiset. Hierholzer finds one cycle in linear time;
# finding them all is a search, and only feasible because the graph is small.
let complete = []
let stack = [{
    at: edges_list[0].to_node,
    used: range(0, len(edges_list)) |> map(|i| i == 0),
    walk: [start_node, edges_list[0].to_node],
}]

while len(stack) > 0 {
    let state = stack[len(stack) - 1]
    stack = slice(stack, 0, len(stack) - 1)

    let remaining = state.used |> count_if(|u| u == false)
    if remaining == 0 {
        # A cycle only counts if it closed.
        if state.at == start_node { complete = push(complete, state.walk) }
    } else {
        for i in range(0, len(edges_list)) {
            if state.used[i] == false and edges_list[i].from_node == state.at {
                let marked = state.used
                marked[i] = true
                stack = push(stack, {
                    at: edges_list[i].to_node,
                    used: marked,
                    walk: push(state.walk, edges_list[i].to_node),
                })
            }
        }
    }
}

# A circular string keeps one character per edge — the trailing k are the leading
# k come round again.
let assembled = complete
    |> map(|walk| range(0, len(walk) - 1) |> map(|i| substr(walk[i], 0, 1)) |> join(""))
    |> unique()
    |> sort()

println("Result:")
for line in assembled { println("  " + line) }
println("Expected: CAGTTCAATTTGGCGTT CAGTTCAATTGGCGTTT CAGTTTCAATTGGCGTT")
println("          CAGTTTGGCGTTCAATT CAGTTGGCGTTCAATTT CAGTTGGCGTTTCAATT")

fn test_grep_assembly_with_repeats() {
    let expected = ["CAGTTCAATTTGGCGTT", "CAGTTCAATTGGCGTTT", "CAGTTTCAATTGGCGTT",
                    "CAGTTTGGCGTTCAATT", "CAGTTGGCGTTCAATTT", "CAGTTGGCGTTTCAATT"]
    assert assembled == sort(expected), "GREP: got " + join(assembled, " ")
    # Each answer is as long as there are reads, and starts with the first one.
    for answer in assembled {
        assert len(answer) == len(patterns),
            "GREP: " + answer + " should be " + str(len(patterns)) + " long"
        assert substr(answer, 0, k + 1) == patterns[0],
            "GREP: " + answer + " must begin with " + patterns[0]
        # And read around the circle, its composition is exactly the input.
        let wrapped = answer + substr(answer, 0, k)
        let composition = range(0, len(answer)) |> map(|i| substr(wrapped, i, k + 1))
        assert sort(composition) == sort(patterns),
            "GREP: " + answer + " does not have the given composition"
    }
    # The point of the problem: more than one genome fits.
    assert len(assembled) > 1, "GREP: repeats should leave the assembly ambiguous"
}

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"
}

OSYM — Isolating Symbols in Alignments

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The best alignment through a given pairing is the best way of reaching it plus the best way of leaving it, so one forward table and one backward table answer all 49 pairs at once — the same trick BA5K uses to find a middle edge.

# Rosalind: OSYM — Isolating Symbols in Alignments
# https://rosalind.info/problems/osym/
#
# Given: Two DNA strings.
# Return: The maximum global alignment score, and the sum over all (j,k) of the
# best score of any alignment that pairs s[j] with t[k].
# Matches score +1, mismatches and gaps -1.

let first_string = "ATAGATA"
let second_string = "ACAGGTA"
let gap = 0 - 1

fn substitution(a, b) { if a == b then 1 else 0 - 1 }

# The question is not what the best alignment is, but how good the best
# alignment *through each particular pairing* would be. Recomputing an alignment
# for every one of the 49 pairs would be quadratic work repeated quadratically.
#
# Instead: the best alignment through (j,k) is the best way of reaching it plus
# the best way of leaving it. One table filled forwards and one filled backwards
# answer every pair at once — the same trick BA5K uses to find a middle edge.
fn forward_table(rows, columns, indel) {
    let table = [range(0, len(columns) + 1) |> map(|j| j * indel)]
    for i in range(1, len(rows) + 1) {
        let line = [i * indel]
        for j in range(1, len(columns) + 1) {
            let diagonal = table[i - 1][j - 1]
                + substitution(substr(rows, i - 1, 1), substr(columns, j - 1, 1))
            let up = table[i - 1][j] + indel
            let left = line[j - 1] + indel
            line = push(line, max([diagonal, up, left]))
        }
        table = push(table, line)
    }
    table
}

let forward = forward_table(first_string, second_string, gap)
# The backward table is the forward one on both strings reversed, so only one
# direction has to be written.
let backward_reversed = forward_table(reverse(first_string), reverse(second_string), gap)

let n = len(first_string)
let m = len(second_string)

fn backward_at(table, i, j, rows, columns) {
    # Score of aligning s[i..) with t[j..), read out of the reversed table.
    table[rows - i][columns - j]
}

let best_score = forward[n][m]

let total = range(1, n + 1) |> flat_map(|j| range(1, m + 1) |> map(|k| {
    let paired = substitution(substr(first_string, j - 1, 1), substr(second_string, k - 1, 1))
    forward[j - 1][k - 1] + paired + backward_at(backward_reversed, j, k, n, m)
})) |> sum()

println("Result:   " + str(best_score))
println("          " + str(total))
println("Expected: 3")
println("          -139")

fn test_osym_isolating_symbols() {
    assert best_score == 3, "OSYM: the best alignment scores " + str(best_score)
    assert total == 0 - 139, "OSYM: the sum is " + str(total)

    # No pairing can beat the unconstrained optimum, and at least one must reach
    # it — the best alignment pairs something.
    let every = range(1, n + 1) |> flat_map(|j| range(1, m + 1) |> map(|k| {
        let paired = substitution(substr(first_string, j - 1, 1), substr(second_string, k - 1, 1))
        forward[j - 1][k - 1] + paired + backward_at(backward_reversed, j, k, n, m)
    }))
    assert max(every) == best_score,
        "OSYM: the best constrained score should equal the unconstrained one"
    for value in every {
        assert value <= best_score, "OSYM: a constrained alignment cannot beat the optimum"
    }
    assert len(every) == n * m, "OSYM: one entry per pair of positions"
}

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"
}

CNTQ — Counting Quartets

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The answer is C(n,4), and the reason matters more than the number: in a fully resolved tree any four taxa are separated into two pairs by some edge, so the shape never enters into it. The assertion enumerates all 15 subsets and confirms it rather than trusting the formula.

# Rosalind: CNTQ — Counting Quartets
# https://rosalind.info/problems/cntq/
#
# Given: n and an unrooted binary tree on n taxa in Newick format.
# Return: The number of quartets consistent with the tree, modulo 1,000,000.

let taxa = ["lobster", "cat", "dog", "caterpillar", "elephant", "mouse"]
# (lobster,(cat,dog),(caterpillar,(elephant,mouse)));
# Each internal edge splits the taxa in two; these are the non-trivial sides.
let clades = [
    ["cat", "dog"],
    ["caterpillar", "elephant", "mouse"],
    ["elephant", "mouse"],
]

# The answer is simply C(n,4), and the reason is worth more than the number: in a
# *fully resolved* binary tree, any four taxa are separated by some internal edge
# into two pairs, so every four-taxon set contributes exactly one consistent
# quartet. Nothing about the tree's shape enters into it.
#
# That stops being true the moment the tree has an unresolved node, which is why
# QRT and CHBP — where splits are partial — are real problems and this one is
# not.
let n = len(taxa)
let answer = choose(n, 4) % 1000000

println("Result:   " + str(answer))
println("Expected: 15")

fn test_cntq_counting_quartets() {
    assert answer == 15, "CNTQ: got " + str(answer)

    # Demonstrated rather than asserted: enumerate all four-taxon subsets and
    # confirm each really is separated into two pairs by some edge.
    fn subsets_of_four(items) {
        range(0, len(items)) |> flat_map(|a|
            range(a + 1, len(items)) |> flat_map(|b|
                range(b + 1, len(items)) |> flat_map(|c|
                    range(c + 1, len(items)) |> map(|d|
                        [items[a], items[b], items[c], items[d]]))))
    }

    let quartets = subsets_of_four(taxa)
    assert len(quartets) == choose(n, 4), "CNTQ: C(6,4) is 15 subsets"

    let consistent = quartets |> count_if(|four|
        (clades |> count_if(|clade|
            (four |> count_if(|t| contains(clade, t))) == 2)) > 0)
    assert consistent == len(quartets),
        "CNTQ: " + str(consistent) + " of " + str(len(quartets)) + " subsets are separated"
    assert consistent == answer, "CNTQ: the count and the formula must agree"
}

QRT — Quartets

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

An 'x' means not scored, not a third state — so a partial character still separates the taxa it did see, and missing data does not make it useless. Quartets are recorded canonically because neither side nor the order within a side carries meaning, and two characters often support the same one.

# Rosalind: QRT — Quartets
# https://rosalind.info/problems/qrt/
#
# Given: A partial character table.
# Return: Every quartet inferable from the splits its characters describe.

let taxa = ["cat", "dog", "elephant", "ostrich", "mouse", "rabbit", "robot"]
let characters = [
    "01xxx00",
    "x11xx00",
    "111x00x",
]

# A partial character is one that could not be scored for every taxon — an 'x'
# means "not known", not "third state". So a character constrains only the taxa
# it actually saw, and the quartets it supports are pairs drawn from its two
# scored groups.
#
# The whole point is that missing data does not make a character useless. It
# still separates the taxa it did score, and those separations are what a tree
# gets built from.
fn group_of(character, symbol, names) {
    range(0, len(names)) |> filter(|i| substr(character, i, 1) == symbol) |> map(|i| names[i])
}

fn pairs_of(items) {
    range(0, len(items)) |> flat_map(|a|
        range(a + 1, len(items)) |> map(|b| [items[a], items[b]]))
}

# A quartet has no orientation — {a,b}|{c,d} is the same as {c,d}|{a,b} — so each
# is recorded in a canonical form and duplicates dropped. Two characters often
# support the same quartet, and counting it twice would overstate the evidence.
fn canonical(left, right) {
    let one = join(sort(left), ",")
    let two = join(sort(right), ",")
    if one < two then one + "|" + two else two + "|" + one
}

let seen = {}
let quartets = []
for character in characters {
    let zeros = group_of(character, "0", taxa)
    let ones = group_of(character, "1", taxa)
    for left in pairs_of(zeros) {
        for right in pairs_of(ones) {
            let key = canonical(left, right)
            if contains(keys(seen), key) == false {
                seen[key] = true
                quartets = push(quartets, { left: left, right: right })
            }
        }
    }
}

let written = quartets
    |> map(|q| "{" + join(q.left, ", ") + "} {" + join(q.right, ", ") + "}")
    |> sort()

println("Result:")
for line in written { println("  " + line) }
println("Expected (any order): {elephant, dog} {rabbit, robot} / {cat, dog} {mouse, rabbit}")
println("                      {mouse, rabbit} {cat, elephant} / {dog, elephant} {mouse, rabbit}")

fn test_qrt_quartets() {
    assert len(quartets) == 4, "QRT: expected 4 quartets, got " + str(len(quartets))
    # Compared as unordered pairs of unordered pairs, since neither side nor the
    # order within a side carries meaning.
    let expected = [
        canonical(["elephant", "dog"], ["rabbit", "robot"]),
        canonical(["cat", "dog"], ["mouse", "rabbit"]),
        canonical(["mouse", "rabbit"], ["cat", "elephant"]),
        canonical(["dog", "elephant"], ["mouse", "rabbit"]),
    ]
    let got = quartets |> map(|q| canonical(q.left, q.right))
    assert sort(got) == sort(expected), "QRT: got " + join(sort(got), " ")
    # Every quartet's four taxa are distinct and really were scored by some
    # character — an 'x' can never appear in one.
    for q in quartets {
        let four = q.left + q.right
        assert len(unique(four)) == 4, "QRT: a quartet must name four different taxa"
        let supporting = characters |> count_if(|c|
            (q.left |> count_if(|t| substr(c, (range(0, len(taxa))
                |> filter(|i| taxa[i] == t))[0], 1) == "0")) == 2
            and (q.right |> count_if(|t| substr(c, (range(0, len(taxa))
                |> filter(|i| taxa[i] == t))[0], 1) == "1")) == 2)
        assert supporting > 0, "QRT: no character supports " + canonical(q.left, q.right)
    }
}

CHBP — Character-Based Phylogeny

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

The inverse of CTBL. It works because a consistent table's splits are laminar — any two nested or disjoint, never crossing — and a laminar family is exactly a tree. Returns a differently-rooted Newick, so the assertion compares induced splits rather than a formatting choice.

# Rosalind: CHBP — Character-Based Phylogeny
# https://rosalind.info/problems/chbp/
#
# Given: Species names and a consistent character table.
# Return: An unrooted binary tree in Newick format modelling the table.

let taxa = ["cat", "dog", "elephant", "mouse", "rabbit", "rat"]
let characters = ["011101", "001101", "001100"]

# The inverse of CTBL: there, a tree produced splits; here the splits have to
# produce the tree. It works because a consistent table's splits are *laminar* —
# any two are nested or disjoint, never crossing — and a laminar family is
# exactly a tree.
#
# An unrooted tree has no first node, so a reference taxon is picked and every
# split is taken from the side away from it. That makes the clades nest, and the
# reference becomes the outermost branch.
let reference = taxa[0]

fn side_away_from(character, names, anchor) {
    let anchor_index = (range(0, len(names)) |> filter(|i| names[i] == anchor))[0]
    let anchor_state = substr(character, anchor_index, 1)
    range(0, len(names))
        |> filter(|i| substr(character, i, 1) != anchor_state)
        |> map(|i| names[i])
}

let clades = characters |> map(|c| side_away_from(c, taxa, reference))

fn is_subset(small, big) { (small |> count_if(|x| contains(big, x) == false)) == 0 }

fn build(members, groups) {
    # The groups strictly inside this one, and of those the maximal — anything
    # contained in another is handled a level deeper.
    let inside = groups |> filter(|g| is_subset(g, members) and len(g) < len(members))
    let maximal = inside |> filter(|g|
        (inside |> count_if(|other| len(other) > len(g) and is_subset(g, other))) == 0)
    let covered = maximal |> flat_map(|g| g)
    let loose = members |> filter(|m| contains(covered, m) == false)
    let parts = (maximal |> map(|g| build(g, groups))) + loose
    if len(parts) == 1 then parts[0] else "(" + join(parts, ",") + ")"
}

let rest = taxa |> filter(|t| t != reference)
let inner = build(rest, clades)
# The reference sits alongside the rest at the unrooted centre, so its branch is
# spliced in rather than wrapped around.
let newick = "(" + reference + "," + substr(inner, 1, len(inner) - 2) + ");"

println("Result:   " + newick)
println("Expected: (dog,(cat,rabbit),(rat,(elephant,mouse)));")
println("Both describe one unrooted tree — they differ only in which branch is")
println("written first, which an unrooted tree does not fix.")

fn test_chbp_character_based_phylogeny() {
    # The real requirement is that the tree induces exactly the table's splits.
    # Newick is not canonical for an unrooted tree, so comparing strings would be
    # comparing a formatting choice.
    fn splits_of(tree_text, names) {
        # Every parenthesised group is a clade; take each as a split.
        let found = []
        let stack = []
        for i in range(0, len(tree_text)) {
            let symbol = substr(tree_text, i, 1)
            if symbol == "(" { stack = push(stack, i) }
            if symbol == ")" {
                let opened = stack[len(stack) - 1]
                stack = slice(stack, 0, len(stack) - 1)
                let inner_text = substr(tree_text, opened + 1, i - opened - 1)
                let members = names |> filter(|t| contains(inner_text, t))
                if len(members) > 1 and len(members) < len(names) {
                    found = push(found, join(sort(members), ","))
                }
            }
        }
        unique(found)
    }

    let mine = splits_of(newick, taxa)
    let published = splits_of("(dog,(cat,rabbit),(rat,(elephant,mouse)));", taxa)

    # A split and its complement are the same split, so compare canonically.
    fn canonical_splits(raw, names) {
        raw |> map(|s| {
            let members = split(s, ",")
            let other = names |> filter(|t| contains(members, t) == false)
            let one = join(sort(members), ",")
            let two = join(sort(other), ",")
            if one < two then one else two
        }) |> unique() |> sort()
    }

    assert canonical_splits(mine, taxa) == canonical_splits(published, taxa),
        "CHBP: splits " + join(canonical_splits(mine, taxa), " | ")
            + " against " + join(canonical_splits(published, taxa), " | ")
    # And those splits are exactly the ones the characters describe.
    let from_characters = canonical_splits(
        clades |> map(|c| join(sort(c), ",")), taxa)
    assert canonical_splits(mine, taxa) == from_characters,
        "CHBP: the tree must induce the table's splits and no others"
}

EUBT — Enumerating Unrooted Binary Trees

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Built by insertion: every unrooted binary tree arises exactly once by splitting an edge and hanging the next taxon off it, which is why the count is (2n-5)!! — fifteen trees at n=5, over two million at n=10. That growth is why nobody enumerates trees to find the best one.

# Rosalind: EUBT — Enumerating Unrooted Binary Trees
# https://rosalind.info/problems/eubt/
#
# Given: n taxa.
# Return: Every unrooted binary tree on those taxa, in Newick format.

let taxa = ["dog", "cat", "mouse", "elephant"]

# Built by insertion: start with the only tree on three taxa — a star — then add
# each remaining taxon by splitting one existing edge in two and hanging it off
# the new node. Every unrooted binary tree arises exactly once this way, which is
# why the count is the double factorial (2n-5)!! and why it explodes: fifteen
# trees at n=5, over two million at n=10.
#
# That growth is the reason nobody enumerates trees to find the best one. It is
# also why the parsimony and distance methods in this pack exist.
let edges_start = [
    { a: "center", b: taxa[0] },
    { a: "center", b: taxa[1] },
    { a: "center", b: taxa[2] },
]

let trees = [edges_start]
for step in range(3, len(taxa)) {
    let leaf = taxa[step]
    let grown = []
    for tree in trees {
        for i in range(0, len(tree)) {
            let fresh = "node" + str(step) + "_" + str(i)
            let split_edge = tree[i]
            let others = range(0, len(tree)) |> filter(|j| j != i) |> map(|j| tree[j])
            grown = push(grown, others + [
                { a: split_edge.a, b: fresh },
                { a: fresh, b: split_edge.b },
                { a: fresh, b: leaf },
            ])
        }
    }
    trees = grown
}

fn neighbours_of(tree, node) {
    (tree |> filter(|e| e.a == node) |> map(|e| e.b))
        + (tree |> filter(|e| e.b == node) |> map(|e| e.a))
}

fn is_taxon(node, names) { contains(names, node) }

fn walk(tree, node, parent, names) {
    if is_taxon(node, names) { return node }
    let children = neighbours_of(tree, node) |> filter(|c| c != parent)
    "(" + (children |> map(|c| walk(tree, c, node, names)) |> join(",")) + ")"
}

# Written from the first taxon, which an unrooted tree does not privilege — it
# is just somewhere to start reading.
fn to_newick(tree, names) {
    let anchor = names[0]
    let hub = neighbours_of(tree, anchor)[0]
    let branches = neighbours_of(tree, hub) |> filter(|c| c != anchor)
    "(" + anchor + "," + (branches |> map(|c| walk(tree, c, hub, names)) |> join(",")) + ");"
}

let written = trees |> map(|t| to_newick(t, taxa))

println("Result:")
for line in written { println("  " + line) }
println("Expected: three trees — (mouse,cat)|(elephant,dog), (elephant,mouse)|(cat,dog),")
println("          (elephant,cat)|(mouse,dog), written from a different anchor.")

fn test_eubt_enumerating_unrooted_binary_trees() {
    # (2n-5)!! trees: 1 x 3 for n = 4.
    assert len(trees) == 3, "EUBT: expected 3 trees, got " + str(len(trees))
    assert len(unique(written)) == 3, "EUBT: the trees must be distinct"

    # A tree on n taxa has 2n-3 edges and n-2 internal nodes.
    for tree in trees {
        assert len(tree) == 2 * len(taxa) - 3,
            "EUBT: expected " + str(2 * len(taxa) - 3) + " edges, got " + str(len(tree))
        for taxon in taxa {
            assert len(neighbours_of(tree, taxon)) == 1, "EUBT: " + taxon + " must be a leaf"
        }
    }

    # The three topologies are exactly the three ways of pairing four taxa. That
    # is the whole content of the answer, and it does not depend on how the
    # Newick is rooted.
    fn pairing_of(tree, names) {
        # The one internal edge separates the taxa into two pairs. Which side the
        # edge happens to be stored from is arbitrary, so the smaller of the two
        # names the topology.
        let internal = tree |> filter(|e|
            is_taxon(e.a, names) == false and is_taxon(e.b, names) == false)
        let side = neighbours_of(tree, internal[0].a) |> filter(|x| is_taxon(x, names))
        let other = names |> filter(|t| contains(side, t) == false)
        let one = join(sort(side), ",")
        let two = join(sort(other), ",")
        if one < two then one else two
    }
    let pairings = trees |> map(|t| pairing_of(t, taxa)) |> sort()
    assert len(unique(pairings)) == 3, "EUBT: each tree must pair the taxa differently"
    # With four taxa the topology is fixed by which one cat is paired with, and
    # all three possibilities appear exactly once.
    assert pairings == sort(["cat,dog", "cat,mouse", "cat,elephant"]),
        "EUBT: got pairings " + join(pairings, " | ")
}

QRTD — Quartet Distance

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Quartets degrade gracefully where splits do not: moving one taxon changes a handful of quartets but can destroy every split at once. Asserted to be zero between a tree and itself, which is the identity any distance must satisfy.

# Rosalind: QRTD — Quartet Distance
# https://rosalind.info/problems/qrtd/
#
# Given: n taxa and two unrooted binary trees.
# Return: The quartet distance between them.

let taxa = ["A", "B", "C", "D", "E"]
# (A,C,((B,D),E));   and   (C,(B,D),(A,E));
# Each internal edge is a split; these are the non-trivial sides.
let first_clades = [["B", "D"], ["B", "D", "E"]]
let second_clades = [["B", "D"], ["A", "E"]]

# Two trees can share every taxon and still disagree about almost everything, so
# comparing them needs a measure. Quartets give one: each set of four taxa is
# separated into two pairs by exactly one edge, and two trees either agree about
# that pairing or they do not.
#
# It beats counting shared splits because it degrades gracefully — moving one
# taxon changes a handful of quartets, where it can destroy every split at once.
fn pairing_in(four, clades) {
    # The clade cutting these four 2-and-2 names the quartet.
    let cutting = clades |> filter(|c| (four |> count_if(|t| contains(c, t))) == 2)
    if len(cutting) == 0 { return "" }
    let side = four |> filter(|t| contains(cutting[0], t))
    let other = four |> filter(|t| contains(cutting[0], t) == false)
    let one = join(sort(side), ",")
    let two = join(sort(other), ",")
    if one < two then one + "|" + two else two + "|" + one
}

fn subsets_of_four(items) {
    range(0, len(items)) |> flat_map(|a|
        range(a + 1, len(items)) |> flat_map(|b|
            range(b + 1, len(items)) |> flat_map(|c|
                range(c + 1, len(items)) |> map(|d|
                    [items[a], items[b], items[c], items[d]]))))
}

let quartets = subsets_of_four(taxa)
let shared = quartets |> count_if(|four| {
    let one = pairing_in(four, first_clades)
    let two = pairing_in(four, second_clades)
    one != "" and one == two
})

# Both trees are fully resolved, so each induces exactly one quartet per subset.
let first_count = quartets |> count_if(|four| pairing_in(four, first_clades) != "")
let second_count = quartets |> count_if(|four| pairing_in(four, second_clades) != "")
let distance = first_count + second_count - 2 * shared

println("Result:   " + str(distance))
println("Expected: 4")
println("(" + str(len(quartets)) + " quartets each, " + str(shared) + " agreeing)")

fn test_qrtd_quartet_distance() {
    assert distance == 4, "QRTD: got " + str(distance)
    # A resolved tree resolves every four-taxon subset, which CNTQ established.
    assert first_count == choose(len(taxa), 4), "QRTD: the first tree resolves all 5"
    assert second_count == choose(len(taxa), 4), "QRTD: and so does the second"
    assert shared == 3, "QRTD: 3 of the 5 quartets agree"

    # A tree against itself is distance zero — the identity any distance must
    # satisfy, and worth checking rather than assuming.
    let self_shared = quartets |> count_if(|four| pairing_in(four, first_clades) != "")
    assert first_count + first_count - 2 * self_shared == 0,
        "QRTD: a tree must be distance 0 from itself"
    # The two disagreeing quartets are the ones involving both B,D and A,E.
    let disagreeing = quartets |> filter(|four|
        pairing_in(four, first_clades) != pairing_in(four, second_clades))
    assert len(disagreeing) == 2, "QRTD: exactly two quartets differ"
}

RNAS — Wobble Bonding and RNA Secondary Structures

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Real RNA pairs U-G nearly as readily as A-U, so a count excluding wobble understates what a molecule can fold into. Memoised with has_key — contains(keys(memo), k) rebuilds the key list on every probe and would make the lookup the bottleneck the memo was meant to remove.

# Rosalind: RNAS — Wobble Bonding and RNA Secondary Structures
# https://rosalind.info/problems/rnas/
#
# Given: An RNA string.
# Return: The number of valid non-crossing matchings of its bonding graph,
# allowing wobble (U-G) pairs, where no bond spans fewer than 4 intervening
# positions.

let rna_string = "AUGCUAGUACGGAGCGAGUCUAGCGAGCGAUGUCGUGAGUACUAUAUAUGCGCAUAAGCCACGU"

# Real RNA does not only pair A-U and C-G. The wobble pair U-G is nearly as
# stable and appears throughout functional RNA, so a structure count that
# excludes it understates what a molecule can fold into.
#
# The other constraint is physical: a hairpin cannot turn in fewer than about
# four bases, so a bond between positions closer than that is impossible however
# well the bases match.
fn can_pair(a, b) {
    (a == "A" and b == "U") or (a == "U" and b == "A")
        or (a == "C" and b == "G") or (a == "G" and b == "C")
        or (a == "U" and b == "G") or (a == "G" and b == "U")
}

# Counting by recursion over intervals: position i either stays unpaired, or
# bonds with some k, which splits the interval into two independent halves. Both
# halves are asked about repeatedly, so without a memo the same interval is
# recomputed exponentially often.
#
# `has_key` is what makes the memo worth having — `contains(keys(memo), k)`
# rebuilds the whole key list on every probe, turning the lookup itself into the
# bottleneck it was meant to remove.
let memo = {}

fn count_matchings(i, j) {
    if j - i < 5 { return 1 }
    let key = str(i) + "," + str(j)
    if has_key(memo, key) { return memo[key] }

    # i unpaired.
    let total = count_matchings(i + 1, j)
    # i bonded to k, which must leave four bases inside the loop.
    for k in range(i + 4, j) {
        if can_pair(substr(rna_string, i, 1), substr(rna_string, k, 1)) {
            total = total + count_matchings(i + 1, k) * count_matchings(k + 1, j)
        }
    }
    memo[key] = total
    total
}

let answer = count_matchings(0, len(rna_string))

println("Result:   " + str(answer))
println("Expected: 284850219977421")

fn test_rnas_wobble_bonding() {
    assert answer == 284850219977421, "RNAS: got " + str(answer)
    # The memo is doing real work: far fewer intervals than the recursion visits.
    assert len(keys(memo)) < len(rna_string) * len(rna_string),
        "RNAS: the memo should hold at most one entry per interval"
    assert len(keys(memo)) > 0, "RNAS: the memo should be used at all"

    # Wobble pairs genuinely matter — U-G is accepted where a strict
    # Watson-Crick rule would reject it.
    assert can_pair("U", "G") and can_pair("G", "U"), "RNAS: wobble pairs bond"
    assert can_pair("A", "G") == false, "RNAS: A-G does not"
    # A string too short to turn a hairpin has exactly one matching: the empty one.
    assert count_matchings(0, 4) == 1, "RNAS: nothing can bond within four bases"
    # The empty matching is always counted, so the total is never zero.
    assert answer > 0, "RNAS: the empty matching always counts"
}

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"
}

RSUB — Identifying Reversing Substitutions

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

A site that mutates and mutates back looks, from the tips alone, as though nothing happened — which is exactly what makes distant relationships hard to recover. Only the internal labels reveal it, which is why the problem supplies them.

# Rosalind: RSUB — Identifying Reversing Substitutions
# https://rosalind.info/problems/rsub/
#
# Given: A rooted binary tree with every node labelled by a string.
# Return: Every reversing substitution — a change that later changes back, with
# no other change in between.

# (((ostrich,cat)rat,mouse)dog,elephant)robot;
let children = {
    "robot": ["dog", "elephant"],
    "dog": ["rat", "mouse"],
    "rat": ["ostrich", "cat"],
}
let sequences = {
    "robot": "AATTG", "dog": "GGGCA", "mouse": "AAGAC", "rat": "GTTGT",
    "cat": "GAGGC", "ostrich": "GTGTC", "elephant": "AATTC",
}
let root = "robot"

# A site that mutates and later mutates back looks, from the tips alone, as
# though nothing ever happened — the ancestor and the descendant agree. Only the
# internal labels reveal it, which is why this problem hands them over rather
# than asking for them.
#
# It matters because reversing substitutions are exactly what makes distant
# relationships hard to recover: the signal erases itself, and two lineages look
# more similar than their history warrants.
fn kids_of(node) { if has_key(children, node) then children[node] else [] }

# Walk down from a node where a change occurred, following only lineages that
# still carry the new character, and report wherever it changes back.
fn reversions_below(start, position, was, became) {
    let found = []
    let frontier = kids_of(start)
    while len(frontier) > 0 {
        let node = frontier[0]
        frontier = slice(frontier, 1, len(frontier))
        let here = substr(sequences[node], position, 1)
        if here == was {
            # Changed back, with nothing else in between.
            found = push(found, node)
        } else {
            # Still carrying the substitution; keep descending. Any third
            # character ends the lineage's relevance.
            if here == became { frontier = frontier + kids_of(node) }
        }
    }
    found
}

let width = len(sequences[root])
let reported = []
let stack = [root]
while len(stack) > 0 {
    let parent = stack[len(stack) - 1]
    stack = slice(stack, 0, len(stack) - 1)
    for child in kids_of(parent) {
        stack = push(stack, child)
        for i in range(0, width) {
            let was = substr(sequences[parent], i, 1)
            let became = substr(sequences[child], i, 1)
            if was != became {
                for reverted in reversions_below(child, i, was, became) {
                    reported = push(reported, child + " " + reverted + " " + str(i + 1)
                        + " " + was + "->" + became + "->" + was)
                }
            }
        }
    }
}

println("Result:")
for line in sort(reported) { println("  " + line) }
println("Expected (any order): dog mouse 1 A->G->A / dog mouse 2 A->G->A")
println("                      rat ostrich 3 G->T->G / rat cat 3 G->T->G / dog rat 3 T->G->T")

fn test_rsub_reversing_substitutions() {
    let expected = ["dog mouse 1 A->G->A", "dog mouse 2 A->G->A",
                    "rat ostrich 3 G->T->G", "rat cat 3 G->T->G", "dog rat 3 T->G->T"]
    assert sort(reported) == sort(expected), "RSUB: got " + join(sort(reported), " / ")

    # Every report must describe a real reversion: the character genuinely
    # differs from the parent and genuinely returns at the named descendant.
    for line in reported {
        let parts = split(line, " ")
        let changed = parts[0]
        let reverted = parts[1]
        let position = int(parts[2]) - 1
        assert substr(sequences[changed], position, 1) != substr(sequences[reverted], position, 1),
            "RSUB: " + line + " does not actually revert"
    }
    # The tips alone hide this: robot and mouse agree at position 1 despite two
    # substitutions on the path between them.
    assert substr(sequences["robot"], 0, 1) == substr(sequences["mouse"], 0, 1),
        "RSUB: the endpoints agree, which is what makes the change invisible"
    assert substr(sequences["dog"], 0, 1) != substr(sequences["robot"], 0, 1),
        "RSUB: yet the intermediate differs from both"
}

REAR — Reversal Distance

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Where BA6C's 2-break distance had a closed form, this has none at size 10 and must be searched: 45 reversals per step, 3.6 million reachable orders, distance reaching 9. The builtin searches from both ends and meets in the middle, so each side only reaches depth 4 or 5.

# Rosalind: REAR — Reversal Distance
# https://rosalind.info/problems/rear/
#
# Given: Up to five pairs of permutations of length 10.
# Return: The reversal distance for each pair.

let pairs = [
    { from_order: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],  to_order: [3, 1, 5, 2, 7, 4, 9, 6, 10, 8] },
    { from_order: [3, 10, 8, 2, 5, 4, 7, 1, 6, 9],  to_order: [5, 2, 3, 1, 7, 4, 10, 8, 6, 9] },
    { from_order: [8, 6, 7, 9, 4, 1, 3, 10, 2, 5],  to_order: [8, 2, 7, 6, 9, 1, 5, 3, 10, 4] },
    { from_order: [3, 9, 10, 4, 1, 8, 6, 7, 5, 2],  to_order: [2, 9, 8, 5, 1, 7, 3, 4, 6, 10] },
    { from_order: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],  to_order: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] },
]

# The 2-break distance in BA6C had a closed form — blocks minus cycles. Reversal
# distance has no such formula at this size, so it has to be searched for, and
# the search is why `reversal_distance` is a builtin rather than written here:
# ten elements admit 45 reversals and 3.6 million reachable orders, and the
# distance reaches 9. A one-sided search to that depth is hopeless; the builtin
# searches from both ends and meets in the middle, so each side only reaches
# depth 4 or 5.
let distances = pairs |> map(|pair| reversal_distance(pair.from_order, pair.to_order))

println("Result:   " + (distances |> map(|d| str(d)) |> join(" ")))
println("Expected: 9 4 5 7 0")

fn test_rear_reversal_distance() {
    assert (distances |> map(|d| str(d)) |> join(" ")) == "9 4 5 7 0",
        "REAR: got " + str(distances)
    # Identical permutations are no distance apart, and the measure is symmetric
    # because a reversal undoes itself.
    assert distances[4] == 0, "REAR: the last pair is already sorted"
    for pair in pairs {
        assert reversal_distance(pair.to_order, pair.from_order)
            == reversal_distance(pair.from_order, pair.to_order),
            "REAR: the distance must be symmetric"
    }
}

SORT — Sorting by Reversals

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

REAR asks how far apart two gene orders are; this asks for the route, which is what actually says how the rearrangement happened. The assertion applies the reversals and checks they produce the target — a plausible list that does not sort is the failure this invites.

# Rosalind: SORT — Sorting by Reversals
# https://rosalind.info/problems/sort/
#
# Given: Two permutations of length 10.
# Return: The reversal distance, and a shortest collection of reversals taking
# the first to the second.

let from_order = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
let to_order = [1, 8, 9, 3, 2, 7, 6, 5, 4, 10]

# REAR asks how far apart two gene orders are; this asks for the route. The
# search is the same — the distance is just the length of what comes back — but
# the sequence is what actually says how the rearrangement happened, and a
# distance without one is unfalsifiable.
#
# Any shortest sequence is correct. Reversals commute in some orders and not
# others, so several routes of the same length usually exist.
let steps = sorting_reversals(from_order, to_order)

println("Result:   " + str(len(steps)))
for step in steps { println("          " + str(step[0]) + " " + str(step[1])) }
println("Expected: 2, then 4 9 and 2 5 — any shortest route is accepted")

fn test_sort_sorting_by_reversals() {
    assert len(steps) == 2, "SORT: expected 2 reversals, got " + str(len(steps))
    assert len(steps) == reversal_distance(from_order, to_order),
        "SORT: the route must be as long as the distance"

    # The route has to work. Applying the reversals in order must produce the
    # target — a plausible-looking list that does not sort is the failure this
    # problem invites.
    let current = from_order
    for step in steps {
        let from_index = step[0] - 1
        let to_index = step[1] - 1
        let head = range(0, from_index) |> map(|i| current[i])
        let middle = range(from_index, to_index + 1)
            |> map(|i| current[to_index - (i - from_index)])
        let tail = range(to_index + 1, len(current)) |> map(|i| current[i])
        current = head + middle + tail
    }
    assert current == to_order,
        "SORT: the reversals give " + str(current) + ", not " + str(to_order)

    # Every interval is a real one, and reverses more than a single element.
    for step in steps {
        assert step[0] >= 1 and step[1] <= len(from_order), "SORT: interval out of range"
        assert step[0] < step[1], "SORT: reversing one element does nothing"
    }
}