Strings

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

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

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

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

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

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

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