Graphs

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

# Rosalind: BA3F — Find an Eulerian Cycle in a Graph
# https://rosalind.info/problems/ba3f/
#
# Given: An Eulerian directed graph, as an adjacency list.
# Return: An Eulerian cycle in the graph.

let adjacency = {
    "0": ["3"],
    "1": ["0"],
    "2": ["1", "6"],
    "3": ["2"],
    "4": ["2"],
    "5": ["4"],
    "6": ["5", "8"],
    "7": ["9"],
    "8": ["7"],
    "9": ["6"],
}

# Hierholzer's algorithm: walk until stuck — which in a balanced graph can only
# happen back where you started — then re-enter at a node with edges left and
# splice the new loop into the walk. Linear in the edges, against the factorial
# cost of searching for the walk directly.
let cycle = eulerian_cycle(adjacency, "6")

println("Result:   " + join(cycle, "->"))
println("Expected: 6->8->7->9->6->5->4->2->1->0->3->2->6")
println("(any Eulerian cycle is accepted — this is one rotation among many)")

fn test_ba3f_eulerian_cycle() {
    assert cycle[0] == cycle[len(cycle) - 1], "BA3F: a cycle must return to its start"

    # The real check is the property, not the published string: every edge used
    # exactly once, and every step a real edge.
    let edge_count = keys(adjacency) |> map(|node| len(adjacency[node])) |> sum()
    assert len(cycle) == edge_count + 1,
        "BA3F: " + str(len(cycle) - 1) + " steps for " + str(edge_count) + " edges"

    let seen = {}
    for i in range(0, len(cycle) - 1) {
        let step = cycle[i] + "->" + cycle[i + 1]
        assert contains(adjacency[cycle[i]], cycle[i + 1]), "BA3F: " + step + " is not an edge"
        assert contains(keys(seen), step) == false, "BA3F: " + step + " is used twice"
        seen[step] = true
    }
    assert len(keys(seen)) == edge_count, "BA3F: not every edge was used"
}

# Rosalind: BA3G — Find an Eulerian Path in a Graph
# https://rosalind.info/problems/ba3g/
#
# Given: A directed graph containing an Eulerian path, as an adjacency list.
# Return: An Eulerian path in the graph.

let adjacency = {
    "0": ["2"],
    "1": ["3"],
    "2": ["1"],
    "3": ["0", "4"],
    "6": ["3", "7"],
    "7": ["8"],
    "8": ["9"],
    "9": ["6"],
}

# A path rather than a cycle, so the walk need not come back. At most one node
# may have an extra edge out — that is where it has to start — and at most one an
# extra edge in, where it has to end. Adding the edge between them makes the
# graph balanced, which turns this into BA3F; the added edge is then cut out
# again, and the walk begins on the far side of the cut.
let path = eulerian_path(adjacency)

println("Result:   " + join(path, "->"))
println("Expected: 6->7->8->9->6->3->0->2->1->3->4")

fn test_ba3g_eulerian_path() {
    # 6 has one more edge out than in, and 4 one more in than out, so the walk is
    # forced to start and end there — that part is not a matter of taste.
    assert path[0] == "6", "BA3G: must start at 6, got " + path[0]
    assert path[len(path) - 1] == "4", "BA3G: must end at 4, got " + path[len(path) - 1]

    let edge_count = keys(adjacency) |> map(|node| len(adjacency[node])) |> sum()
    assert len(path) == edge_count + 1,
        "BA3G: " + str(len(path) - 1) + " steps for " + str(edge_count) + " edges"

    let seen = {}
    for i in range(0, len(path) - 1) {
        let step = path[i] + "->" + path[i + 1]
        assert contains(adjacency[path[i]], path[i + 1]), "BA3G: " + step + " is not an edge"
        assert contains(keys(seen), step) == false, "BA3G: " + step + " is used twice"
        seen[step] = true
    }
    assert len(keys(seen)) == edge_count, "BA3G: not every edge was used"
}

# Rosalind: BA5N — Find a Topological Ordering of a DAG
# https://rosalind.info/problems/ba5n/
#
# Given: The adjacency list of a directed acyclic graph.
# Return: A topological ordering of its vertices.

let adjacency = {
    "1": ["2"],
    "2": ["3"],
    "4": ["2"],
    "5": ["3"],
}

# An order in which every edge points forwards. This is what makes BA5B and BA5D
# possible: a longest path can be found in one sweep only if each node is reached
# after everything that leads into it.
#
# Kahn's algorithm: repeatedly take a node nothing points at, remove it, and see
# what that frees. The initial set is sorted so the answer does not depend on
# hash iteration order; freed vertices then join the back of the queue, which is
# what makes this first-in-first-out rather than a re-sort each round.
#
# Any ordering with every edge pointing forwards is correct — several exist here,
# and the assertion checks that property rather than only the printed string.
let vertices = sort(unique(keys(adjacency) + (keys(adjacency) |> flat_map(|k| adjacency[k]))))

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

let ready = vertices |> filter(|node| incoming[node] == 0) |> sort()
let ordering = []
while len(ready) > 0 {
    let node = ready[0]
    ready = slice(ready, 1, len(ready))
    ordering = push(ordering, node)
    if contains(keys(adjacency), node) {
        for target in adjacency[node] {
            incoming[target] = incoming[target] - 1
            if incoming[target] == 0 { ready = push(ready, target) }
        }
    }
}

println("Result:   " + join(ordering, ", "))
println("Expected: 1, 4, 5, 2, 3")

fn test_ba5n_topological_ordering() {
    assert join(ordering, ", ") == "1, 4, 5, 2, 3", "BA5N: got " + join(ordering, ", ")
    # Every node appears exactly once...
    assert sort(ordering) == vertices, "BA5N: the ordering must list every node once"
    # ...and every edge points forwards, which is the whole definition.
    let position = {}
    for i in range(0, len(ordering)) { position[ordering[i]] = i }
    for node in keys(adjacency) {
        for target in adjacency[node] {
            assert position[node] < position[target],
                "BA5N: " + node + " -> " + target + " points backwards"
        }
    }
}

# Rosalind: BA3M — Generate All Maximal Non-Branching Paths in a Graph
# https://rosalind.info/problems/ba3m/
#
# Given: The adjacency list of a directed graph.
# Return: Every maximal non-branching path.

let adjacency = {
    "1": ["2"], "2": ["3"], "3": ["4", "5"], "6": ["7"], "7": ["6"],
}

# A run of nodes with one edge in and one out carries no choice — nothing about
# the graph is lost by collapsing it to a single path. Everywhere else there is a
# genuine branch, and that is where an assembly has to stop and admit it does not
# know which way the genome went. This is the operation that turns a de Bruijn
# graph into contigs.
let vertices = sort(unique(keys(adjacency) + (keys(adjacency) |> flat_map(|k| adjacency[k]))))

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(adjacency) {
    for target in adjacency[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
}

let paths = []
let used = {}
for node in vertices {
    if is_one_in_one_out(node, adjacency, in_degree) == false {
        for target in out_edges(node, adjacency) {
            let walk = [node, target]
            used[node + "->" + target] = true
            let at = target
            while is_one_in_one_out(at, adjacency, in_degree) {
                let onward = out_edges(at, adjacency)[0]
                used[at + "->" + onward] = true
                walk = push(walk, onward)
                at = onward
            }
            paths = push(paths, walk)
        }
    }
}

# Isolated cycles: every node 1-in-1-out, so no starting point was ever found.
for node in vertices {
    if is_one_in_one_out(node, adjacency, in_degree) {
        let first_edge = node + "->" + out_edges(node, adjacency)[0]
        if contains(keys(used), first_edge) == false {
            let walk = [node]
            let at = node
            let going = true
            while going {
                let onward = out_edges(at, adjacency)[0]
                used[at + "->" + onward] = true
                walk = push(walk, onward)
                at = onward
                if at == node { going = false }
            }
            paths = push(paths, walk)
        }
    }
}

let listed = sort(paths |> map(|p| join(p, " -> ")))

println("Result:")
for line in listed { println("  " + line) }
println("Expected: 1 -> 2 -> 3 / 3 -> 4 / 3 -> 5 / 7 -> 6 -> 7")
println("(the cycle is printed from 6 rather than 7 — a cycle has no first node)")

fn test_ba3m_maximal_non_branching_paths() {
    # The three linear paths are pinned down exactly; the cycle is only pinned
    # up to where it starts, because 6 -> 7 -> 6 and 7 -> 6 -> 7 are the same
    # cycle traversed from different nodes.
    for wanted in ["1 -> 2 -> 3", "3 -> 4", "3 -> 5"] {
        assert contains(listed, wanted), "BA3M: missing " + wanted
    }
    assert len(listed) == 4, "BA3M: expected 4 paths, got " + str(len(listed))
    let cycles = paths |> filter(|p| p[0] == p[len(p) - 1])
    assert len(cycles) == 1, "BA3M: exactly one isolated cycle"
    assert sort(unique(cycles[0])) == ["6", "7"], "BA3M: the cycle runs through 6 and 7"
    # Every edge belongs to exactly one path — the paths partition the graph,
    # which is what makes them a lossless summary of it.
    let edge_count = keys(adjacency) |> map(|node| len(adjacency[node])) |> sum()
    let covered = paths |> map(|p| len(p) - 1) |> sum()
    assert covered == edge_count,
        "BA3M: paths cover " + str(covered) + " edges of " + str(edge_count)
    # The 6-7 cycle has no non-branching start, so it is only found by the second
    # pass — without which it would be silently dropped.
    assert len(cycles) == 1, "BA3M: the isolated cycle must be reported"
}