Combinatorics

15 problems from Rosalind — Bioinformatics Stronghold. Press Run on any block to execute it in your browser.

# 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)
}

# 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)
}

# 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)
}

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

# 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)
}

# 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)
}

# 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)
}

# 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)
}

# 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)
}

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

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

# 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]
}

# 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)
}

# 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)
}

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