Assembly

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

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

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

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

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

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