Assembly

9 problems from Rosalind — Bioinformatics Textbook Track. Press Run on any block to execute it in your browser.

# Rosalind: BA3B — Reconstruct a String from its Genome Path
# https://rosalind.info/problems/ba3b/
#
# Given: A sequence of k-mers where each overlaps the next by k-1 symbols.
# Return: The string they spell.

let path = ["ACCGA", "CCGAA", "CGAAG", "GAAGC", "AAGCT"]

# The path is already in order, so the string is the first k-mer plus the last
# symbol of each one after it. No overlap has to be searched for — that is what
# BA3C and BA3H are for.
let text = range(1, len(path))
  |> reduce(|acc, i| acc + substr(path[i], len(path[i]) - 1, 1), path[0])

println("Result:   " + text)
println("Expected: ACCGAAGCT")

fn test_ba3b_genome_path() {
    assert text == "ACCGAAGCT", "BA3B: got " + text
    assert len(text) == len(path[0]) + len(path) - 1, "BA3B: length should be k + n - 1"
}

# Rosalind: BA3C — Construct the Overlap Graph of a Collection of k-mers
# https://rosalind.info/problems/ba3c/
#
# Given: A collection of k-mers.
# Return: The overlap graph as an adjacency list.

let patterns = ["ATGCG", "GCATG", "CATGC", "AGGCA", "GGCAT"]

# An edge runs from one k-mer to another when the suffix of the first is the
# prefix of the second — the relation assembly walks. A k-mer is not joined to
# itself unless it genuinely overlaps itself.
fn suffix_of(pattern) { substr(pattern, 1, len(pattern) - 1) }
fn prefix_of(pattern) { substr(pattern, 0, len(pattern) - 1) }

let overlaps = []
for source in patterns {
    for target in patterns {
        if source != target and suffix_of(source) == prefix_of(target) {
            overlaps = push(overlaps, source + " -> " + target)
        }
    }
}

let result = sort(overlaps)

println("Result:")
for edge in result {
    println("  " + edge)
}
println("Expected: AGGCA -> GGCAT / CATGC -> ATGCG / GCATG -> CATGC / GGCAT -> GCATG")

fn test_ba3c_overlap_graph() {
    assert len(result) == 4, "BA3C: expected 4 overlaps, got " + str(len(result))
    for expected in ["AGGCA -> GGCAT", "CATGC -> ATGCG", "GCATG -> CATGC", "GGCAT -> GCATG"] {
        assert contains(result, expected), "BA3C: missing " + expected
    }
}

# Rosalind: BA3D — Construct the De Bruijn Graph of a String
# https://rosalind.info/problems/ba3d/
#
# Given: An integer k and a string Text.
# Return: DeBruijn_k(Text) as an adjacency list.

let k = 4
let text = "AAGATTCTCTAC"

# Nodes are (k-1)-mers and edges are the k-mers: each k-mer joins its own prefix
# to its own suffix. That is the inversion that makes assembly an Eulerian path
# problem rather than a Hamiltonian one — edges are what must be used, not nodes.
let adjacency = {}
for i in range(0, len(text) - k + 1) {
    let mer = substr(text, i, k)
    let source = substr(mer, 0, k - 1)
    let target = substr(mer, 1, k - 1)
    if contains(keys(adjacency), source) {
        adjacency[source] = push(adjacency[source], target)
    } else {
        adjacency[source] = [target]
    }
}

# Repeated edges are kept — a k-mer occurring twice is two edges, which is what
# lets the assembly traverse a repeat the right number of times.
let result = sort(keys(adjacency)) |> map(|node| node + " -> " + (sort(adjacency[node]) |> join(",")))

println("Result:")
for line in result {
    println("  " + line)
}
println("Expected: AAG -> AGA / TCT -> CTA,CTC / ... (8 lines)")

fn test_ba3d_de_bruijn_of_a_string() {
    assert len(result) == 8, "BA3D: expected 8 nodes, got " + str(len(result))
    assert contains(result, "TCT -> CTA,CTC"), "BA3D: TCT should branch to both"
    assert contains(result, "AAG -> AGA"), "BA3D: missing AAG -> AGA"
    # Every k-mer of the text is one edge.
    let edge_count = result |> map(|l| len(split(substr(l, index_of(l, "-> ") + 3, len(l)), ","))) |> sum()
    assert edge_count == len(text) - k + 1, "BA3D: edges should equal the number of k-mers"
}

# Rosalind: BA3E — Construct the De Bruijn Graph of a Collection of k-mers
# https://rosalind.info/problems/ba3e/
#
# Given: A collection of k-mers.
# Return: The de Bruijn graph as an adjacency list.

let patterns = ["GAGG", "CAGG", "GGGG", "GGGA", "CAGG", "AGGG", "GGAG"]

# Same construction as BA3D, but the k-mers arrive as a bag rather than being
# read out of a string — which is the realistic case, since reads are what a
# sequencer produces. CAGG appears twice in the input and must appear twice in
# the graph: duplicate reads are evidence of a repeat, not noise to discard.
let adjacency = {}
for mer in patterns {
    let source = substr(mer, 0, len(mer) - 1)
    let target = substr(mer, 1, len(mer) - 1)
    if contains(keys(adjacency), source) {
        adjacency[source] = push(adjacency[source], target)
    } else {
        adjacency[source] = [target]
    }
}

let result = sort(keys(adjacency)) |> map(|node| node + " -> " + (sort(adjacency[node]) |> join(",")))

println("Result:")
for line in result {
    println("  " + line)
}
println("Expected: AGG -> GGG / CAG -> AGG,AGG / GAG -> AGG / GGA -> GAG / GGG -> GGA,GGG")

fn test_ba3e_de_bruijn_of_a_collection() {
    assert len(result) == 5, "BA3E: expected 5 nodes, got " + str(len(result))
    assert contains(result, "CAG -> AGG,AGG"), "BA3E: the repeated k-mer must stay repeated"
    assert contains(result, "GGG -> GGA,GGG"), "BA3E: missing GGG's two edges"
}

# Rosalind: BA3H — Reconstruct a String from its k-mer Composition
# https://rosalind.info/problems/ba3h/
#
# Given: An integer k, followed by a list of k-mers Patterns.
# Return: A string Text whose k-mer composition is Patterns.

let k = 4
let patterns = ["CTTA", "ACCA", "TACC", "GGCT", "GCTT", "TTAC"]

# Genome assembly, in miniature. Each k-mer becomes an *edge* from its prefix to
# its suffix — not a node — because then using every read exactly once is
# precisely an Eulerian path, which BA3G already solves in linear time. Making
# reads the nodes instead gives a Hamiltonian path, which is NP-hard: the same
# data, a different graph, and the difference between tractable and not.
let de_bruijn = {}
for pattern in patterns {
    let prefix = substr(pattern, 0, k - 1)
    let suffix = substr(pattern, 1, k - 1)
    if contains(keys(de_bruijn), prefix) {
        de_bruijn[prefix] = push(de_bruijn[prefix], suffix)
    } else {
        de_bruijn[prefix] = [suffix]
    }
}

let path = eulerian_path(de_bruijn)

# Spell the path: the first node in full, then one new character per step.
let text = path[0] + (range(1, len(path))
    |> map(|i| substr(path[i], k - 2, 1))
    |> join(""))

println("Result:   " + text)
println("Expected: GGCTTACCA")

fn test_ba3h_string_reconstruction() {
    assert text == "GGCTTACCA", "BA3H: got " + text
    assert len(text) == len(patterns) + k - 1, "BA3H: n reads of length k spell n + k - 1"
    # The real requirement: the answer's k-mer composition is the input, as a
    # multiset. Sorting both is enough since every read is used once.
    let composition = range(0, len(text) - k + 1) |> map(|i| substr(text, i, k))
    assert sort(composition) == sort(patterns),
        "BA3H: composition " + str(sort(composition)) + " != " + str(sort(patterns))
}

# Rosalind: BA3I — Find a k-Universal Circular String
# https://rosalind.info/problems/ba3i/
#
# Given: An integer k.
# Return: A k-universal circular binary string — one containing every binary
# k-mer exactly once when read around the circle.

let k = 4

# Every binary k-mer, in order: 0000, 0001, ... 1111.
fn binary_kmers(width) {
    let patterns = [""]
    for _ in range(0, width) {
        patterns = patterns |> flat_map(|prefix| [prefix + "0", prefix + "1"])
    }
    patterns
}

let patterns = binary_kmers(k)

# The same de Bruijn construction as BA3H, but here every (k-1)-mer has exactly
# two edges in and two out, so the graph is balanced and the walk closes into a
# cycle. A circular string of length 2^k containing all 2^k k-mers is only
# possible because each one overlaps the next by k-1 — the cycle is what makes
# that packing exist at all.
let de_bruijn = {}
for pattern in patterns {
    let prefix = substr(pattern, 0, k - 1)
    let suffix = substr(pattern, 1, k - 1)
    if contains(keys(de_bruijn), prefix) {
        de_bruijn[prefix] = push(de_bruijn[prefix], suffix)
    } else {
        de_bruijn[prefix] = [suffix]
    }
}

let cycle = eulerian_cycle(de_bruijn, patterns[0] |> substr(0, k - 1))

# Spelling a *circular* string drops the last k-1 characters, because they are
# the first k-1 come round again.
let spelled = cycle[0] + (range(1, len(cycle))
    |> map(|i| substr(cycle[i], k - 2, 1))
    |> join(""))
let text = substr(spelled, 0, len(spelled) - (k - 1))

println("Result:   " + text)
println("Expected: 0000110010111101  (any k-universal string is accepted)")

fn test_ba3i_k_universal_circular_string() {
    assert len(text) == pow(2, k), "BA3I: a k-universal binary string has length 2^k"
    # Read around the circle: every binary k-mer exactly once.
    let wrapped = text + substr(text, 0, k - 1)
    let found = range(0, len(text)) |> map(|i| substr(wrapped, i, k))
    assert len(unique(found)) == pow(2, k), "BA3I: a k-mer appears twice"
    assert sort(found) == sort(patterns), "BA3I: not every k-mer appears"
}

# Rosalind: BA3J — Reconstruct a String from its Paired Composition
# https://rosalind.info/problems/ba3j/
#
# Given: Integers k and d, followed by a collection of (k,d)-mers.
# Return: A string whose (k,d)-mer composition is the given collection.

let k = 4
let d = 2
let pairs = [
    "GAGA|TTGA", "TCGT|GATG", "CGTG|ATGT", "TGGT|TGAG", "GTGA|TGTT",
    "GTGG|GTGA", "TGAG|GTTG", "GGTC|GAGA", "GTCG|AGAT",
]

# Paired reads carry information plain k-mers do not: two short reads a known
# distance apart pin down a repeat that either one alone would be ambiguous
# inside. The graph is built the same way as BA3H, on pairs rather than single
# k-mers — the prefix of a pair is the prefix of both halves.
fn halves(pair) { split(pair, "|") }

let de_bruijn = {}
for pair in pairs {
    let parts = halves(pair)
    let prefix = substr(parts[0], 0, k - 1) + "|" + substr(parts[1], 0, k - 1)
    let suffix = substr(parts[0], 1, k - 1) + "|" + substr(parts[1], 1, k - 1)
    if contains(keys(de_bruijn), prefix) {
        de_bruijn[prefix] = push(de_bruijn[prefix], suffix)
    } else {
        de_bruijn[prefix] = [suffix]
    }
}

let path = eulerian_path(de_bruijn)

# Spell each half separately, then overlap them. The two spellings must agree
# wherever they overlap — that agreement is the extra constraint the pairing
# buys, and checking it is what makes a wrong assembly detectable.
fn spell(nodes, which, width) {
    let first = halves(nodes[0])[which]
    first + (range(1, len(nodes)) |> map(|i| substr(halves(nodes[i])[which], width - 2, 1)) |> join(""))
}

let first_spelled = spell(path, 0, k)
let second_spelled = spell(path, 1, k)
let gap = k + d

let overlap_disagrees = range(gap, len(first_spelled))
    |> filter(|i| substr(first_spelled, i, 1) != substr(second_spelled, i - gap, 1))

let text = first_spelled + substr(second_spelled, len(second_spelled) - gap, gap)

println("Result:   " + text)
println("Expected: GTGGTCGTGAGATGTTGA")

fn test_ba3j_paired_reconstruction() {
    assert text == "GTGGTCGTGAGATGTTGA", "BA3J: got " + text
    assert len(overlap_disagrees) == 0,
        "BA3J: the two spellings disagree at " + str(overlap_disagrees)
    assert len(text) == len(pairs) + 2 * k + d - 1,
        "BA3J: n pairs spell n + 2k + d - 1 characters"
    # And the answer's own paired composition is the input.
    let composition = range(0, len(text) - (2 * k + d) + 1)
        |> map(|i| substr(text, i, k) + "|" + substr(text, i + k + d, k))
    assert sort(composition) == sort(pairs),
        "BA3J: composition does not match the reads"
}

# Rosalind: BA3K — Generate Contigs from a Collection of Reads
# https://rosalind.info/problems/ba3k/
#
# Given: A collection of k-mers Patterns.
# Return: The contigs from the de Bruijn graph of Patterns.

let patterns = ["ATG", "ATG", "TGT", "TGG", "CAT", "GGA", "GAT", "AGA"]
let k = 3

# What assembly actually produces. BA3H asked for *the* genome, which needs an
# Eulerian path to exist and be unique — real data gives neither. Wherever the
# graph branches, the reads genuinely do not say which way the genome went, so
# the honest output is the unambiguous stretches and no more. Those are the
# contigs, and their lengths are the standard measure of how good an assembly is.
let de_bruijn = {}
for pattern in patterns {
    let prefix = substr(pattern, 0, k - 1)
    let suffix = substr(pattern, 1, k - 1)
    if contains(keys(de_bruijn), prefix) {
        de_bruijn[prefix] = push(de_bruijn[prefix], suffix)
    } else {
        de_bruijn[prefix] = [suffix]
    }
}

let vertices = sort(unique(keys(de_bruijn) + (keys(de_bruijn) |> flat_map(|n| de_bruijn[n]))))
fn out_edges(node, graph) { if contains(keys(graph), node) then graph[node] else [] }

let in_degree = {}
for node in vertices { in_degree[node] = 0 }
for node in keys(de_bruijn) {
    for target in de_bruijn[node] { in_degree[target] = in_degree[target] + 1 }
}

fn is_one_in_one_out(node, graph, degrees) {
    degrees[node] == 1 and len(out_edges(node, graph)) == 1
}

# Maximal non-branching paths, exactly as in BA3M.
let paths = []
for node in vertices {
    if is_one_in_one_out(node, de_bruijn, in_degree) == false {
        for target in out_edges(node, de_bruijn) {
            let walk = [node, target]
            let at = target
            while is_one_in_one_out(at, de_bruijn, in_degree) {
                let onward = out_edges(at, de_bruijn)[0]
                walk = push(walk, onward)
                at = onward
            }
            paths = push(paths, walk)
        }
    }
}

# Spell each path, as in BA3H.
let contigs = sort(paths |> map(|walk|
    walk[0] + (range(1, len(walk)) |> map(|i| substr(walk[i], k - 2, 1)) |> join(""))))

println("Result:   " + join(contigs, " "))
println("Expected: AGA ATG ATG CAT GAT TGGA TGT")

fn test_ba3k_contig_generation() {
    assert join(contigs, " ") == "AGA ATG ATG CAT GAT TGGA TGT",
        "BA3K: got " + join(contigs, " ")
    # ATG appears twice, because two separate branches spell it — contigs are
    # listed per path, not deduplicated.
    assert (contigs |> count_if(|c| c == "ATG")) == 2, "BA3K: ATG is spelled by two paths"
    # TGGA is the only contig longer than k, which is exactly where the graph did
    # not branch. Everything else stops immediately.
    assert (contigs |> count_if(|c| len(c) > k)) == 1, "BA3K: only one contig extends past k"
    assert contains(contigs, "TGGA"), "BA3K: TGGA is the unambiguous stretch"
    # Every contig is spellable from the reads it came from.
    for contig in contigs {
        let pieces = range(0, len(contig) - k + 1) |> map(|i| substr(contig, i, k))
        for piece in pieces {
            assert contains(patterns, piece), "BA3K: " + contig + " uses a read that is not there"
        }
    }
}

# Rosalind: BA3L — Construct a String Spelled by a Gapped Genome Path
# https://rosalind.info/problems/ba3l/
#
# Given: A sequence of (k,d)-mers already in path order.
# Return: The string they spell, if one exists.

let k = 4
let d = 2
let path = [
    "GACC|GCGC", "ACCG|CGCC", "CCGA|GCCG", "CGAG|CCGG", "GAGC|CGGA",
]

# BA3J had to find the path first; here it is given, so what remains is the
# spelling — and the check that the path is consistent at all. The two halves
# overlap by k+d characters once laid out, and if they disagree anywhere then no
# string has this paired composition, however valid the path looked in the graph.
fn halves(pair) { split(pair, "|") }

fn spell(nodes, which, width) {
    let head = halves(nodes[0])[which]
    head + (range(1, len(nodes)) |> map(|i| substr(halves(nodes[i])[which], width - 1, 1)) |> join(""))
}

let first_spelled = spell(path, 0, k)
let second_spelled = spell(path, 1, k)
let gap = k + d

let disagreements = range(gap, len(first_spelled))
    |> filter(|i| substr(first_spelled, i, 1) != substr(second_spelled, i - gap, 1))

let text = first_spelled + substr(second_spelled, len(second_spelled) - gap, gap)

println("Result:   " + text)
println("Expected: GACCGAGCGCCGGA")

fn test_ba3l_gapped_genome_path() {
    assert text == "GACCGAGCGCCGGA", "BA3L: got " + text
    # The overlap check is the substance of this problem, not a formality: a path
    # whose halves disagree spells nothing at all.
    assert len(disagreements) == 0,
        "BA3L: the two spellings disagree at " + str(disagreements)
    assert len(text) == len(path) + 2 * k + d - 1,
        "BA3L: n pairs spell n + 2k + d - 1 characters"
    # And the answer's own paired composition is the path it came from.
    let composition = range(0, len(text) - (2 * k + d) + 1)
        |> map(|i| substr(text, i, k) + "|" + substr(text, i + k + d, k))
    assert composition == path, "BA3L: the composition does not match the given path"
}