Rearrangements

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

# Rosalind: BA6A — Implement GreedySorting to Sort a Permutation by Reversals
# https://rosalind.info/problems/ba6a/
#
# Given: A signed permutation P.
# Return: Every permutation produced by GreedySorting, ending at the identity.

let start = [-3, 4, 1, 5, -2]

# Chromosomes rearrange by reversal — a segment is cut out, flipped and put back,
# which is why the signs matter: a reversed gene reads on the other strand. The
# number of reversals separating two genomes is a measure of how long ago they
# diverged.
#
# Greedy sorting fixes position 1 first, then 2, and never revisits either. It
# takes at most 2n reversals, which is not the minimum — finding that is much
# harder — but it terminates and it bounds the true distance.
fn reverse_segment(items, from_index, to_index) {
    let head = range(0, from_index) |> map(|i| items[i])
    let middle = range(from_index, to_index + 1) |> map(|i| 0 - items[to_index - (i - from_index)])
    let tail = range(to_index + 1, len(items)) |> map(|i| items[i])
    head + middle + tail
}

let permutation = start
let steps = []
for k in range(0, len(start)) {
    if permutation[k] != k + 1 {
        # Where the value that belongs here currently sits, either sign.
        let at = (range(k, len(permutation)) |> filter(|i| abs(permutation[i]) == k + 1))[0]
        permutation = reverse_segment(permutation, k, at)
        steps = push(steps, permutation)
        # The reversal may have left it negative, which costs one more flip.
        if permutation[k] == 0 - (k + 1) {
            permutation = reverse_segment(permutation, k, k)
            steps = push(steps, permutation)
        }
    }
}

fn written(items) {
    "(" + (items |> map(|v| if v > 0 then "+" + str(v) else str(v)) |> join(" ")) + ")"
}

println("Result:")
for step in steps { println("  " + written(step)) }
println("Expected: (-1 -4 +3 +5 -2) (+1 -4 +3 +5 -2) (+1 +2 -5 -3 +4)")
println("          (+1 +2 +3 +5 +4) (+1 +2 +3 -4 -5) (+1 +2 +3 +4 -5) (+1 +2 +3 +4 +5)")

fn test_ba6a_greedy_sorting() {
    let expected = ["(-1 -4 +3 +5 -2)", "(+1 -4 +3 +5 -2)", "(+1 +2 -5 -3 +4)",
                    "(+1 +2 +3 +5 +4)", "(+1 +2 +3 -4 -5)", "(+1 +2 +3 +4 -5)",
                    "(+1 +2 +3 +4 +5)"]
    assert (steps |> map(|s| written(s))) == expected,
        "BA6A: got " + join(steps |> map(|s| written(s)), " ")
    # It really ends at the identity, and takes at most 2n reversals.
    assert steps[len(steps) - 1] == range(1, len(start) + 1), "BA6A: must end sorted"
    assert len(steps) <= 2 * len(start), "BA6A: greedy sorting is bounded by 2n reversals"
    # Every step is a genuine reversal of the one before: same values, ignoring
    # sign and order.
    let running = start
    for step in steps {
        assert sort(step |> map(|v| abs(v))) == sort(running |> map(|v| abs(v))),
            "BA6A: a reversal cannot change which values are present"
        running = step
    }
}

# Rosalind: BA6B — Compute the Number of Breakpoints in a Permutation
# https://rosalind.info/problems/ba6b/
#
# Given: A signed permutation P.
# Return: The number of breakpoints in P.

let permutation = [3, 4, 5, -12, -8, -7, -6, 1, 2, 10, 9, -11, 13, 14]

# A breakpoint is any adjacent pair that is not consecutive — anywhere the
# permutation is already in order, no reversal need ever touch. Since one
# reversal can remove at most two breakpoints, the breakpoint count divided by
# two is a lower bound on the true reversal distance, which is what makes this
# worth computing at all: BA6A gives an upper bound, this gives a lower one.
#
# The permutation is bracketed by 0 and n+1 so the ends count too — a chromosome
# whose first gene is not gene 1 has a breakpoint there.
let extended = [0] + permutation + [len(permutation) + 1]

let breakpoints = range(1, len(extended))
    |> count_if(|i| extended[i] - extended[i - 1] != 1)

println("Result:   " + str(breakpoints))
println("Expected: 8")

fn test_ba6b_breakpoints() {
    assert breakpoints == 8, "BA6B: got " + str(breakpoints)
    # The identity permutation has none, which is the only permutation that does.
    let identity = [0] + range(1, 6) + [6]
    assert (range(1, len(identity)) |> count_if(|i| identity[i] - identity[i - 1] != 1)) == 0,
        "BA6B: the identity has no breakpoints"
    # Worth being careful here: a *fully reversed* permutation is not the worst
    # case. (-5 -4 -3 -2 -1) still steps by one at every position, so it has only
    # the two breakpoints at its ends — and one reversal sorts it, which the
    # bound below then predicts exactly.
    let flipped = [0, -5, -4, -3, -2, -1, 6]
    assert (range(1, len(flipped)) |> count_if(|i| flipped[i] - flipped[i - 1] != 1)) == 2,
        "BA6B: a full reversal leaves only the two end breakpoints"
    # Interleaving is what actually breaks every adjacency.
    let shuffled = [0, 2, 4, 1, 3, 5]
    assert (range(1, len(shuffled)) |> count_if(|i| shuffled[i] - shuffled[i - 1] != 1)) == 5,
        "BA6B: interleaving breaks every one of the five adjacencies"
    # The lower bound this implies on reversal distance.
    assert breakpoints / 2 == 4, "BA6B: at least 4 reversals are needed"
}

# Rosalind: BA6C — Compute the 2-Break Distance Between a Pair of Genomes
# https://rosalind.info/problems/ba6c/
#
# Given: Two genomes with circular chromosomes on the same synteny blocks.
# Return: The 2-break distance between them.

let genome_p = [[1, 2, 3, 4, 5, 6]]
let genome_q = [[1, -3, -6, -5], [2, -4]]

# There is a closed form, which is unusual for a rearrangement distance and is
# the reason 2-breaks are studied at all: the distance is the number of synteny
# blocks minus the number of cycles in the two genomes' graphs superimposed.
#
# Superimposing them, every node has one edge from each genome, so the graph
# splits into alternating cycles. A 2-break can increase the cycle count by at
# most one, and the genomes are identical exactly when every cycle is trivial —
# so blocks minus cycles is both a lower bound and achievable.
fn coloured_edges(genome) {
    genome |> flat_map(|chromosome| {
        let nodes = chromosome |> flat_map(|block|
            if block > 0 then [2 * block - 1, 2 * block] else [0 - 2 * block, 0 - 2 * block - 1])
        range(1, len(chromosome) + 1) |> map(|j| {
            let to_index = if 2 * j == len(nodes) then 0 else 2 * j
            { a: nodes[2 * j - 1], b: nodes[to_index] }
        })
    })
}

let edges_p = coloured_edges(genome_p)
let edges_q = coloured_edges(genome_q)
let blocks = genome_p |> map(|c| len(c)) |> sum()

# Adjacency over both edge sets at once.
let links = {}
for edge in edges_p + edges_q {
    for pair in [{ x: edge.a, y: edge.b }, { x: edge.b, y: edge.a }] {
        if contains(keys(links), str(pair.x)) {
            links[str(pair.x)] = push(links[str(pair.x)], pair.y)
        } else {
            links[str(pair.x)] = [pair.y]
        }
    }
}

let seen = {}
let cycles = 0
for node in keys(links) {
    if contains(keys(seen), node) == false {
        cycles = cycles + 1
        let frontier = [node]
        while len(frontier) > 0 {
            let at = frontier[0]
            frontier = slice(frontier, 1, len(frontier))
            if contains(keys(seen), at) == false {
                seen[at] = true
                for neighbour in links[at] {
                    if contains(keys(seen), str(neighbour)) == false {
                        frontier = push(frontier, str(neighbour))
                    }
                }
            }
        }
    }
}

let distance = blocks - cycles

println("Result:   " + str(distance))
println("Expected: 3")
println("(" + str(blocks) + " blocks minus " + str(cycles) + " cycles)")

fn test_ba6c_two_break_distance() {
    assert distance == 3, "BA6C: got " + str(distance)
    assert blocks == 6, "BA6C: six synteny blocks"
    assert cycles == 3, "BA6C: the superimposed graph has three cycles"
    # A genome against itself is distance zero, which is the identity the formula
    # has to satisfy: every cycle is trivial, so cycles equal blocks.
    let self_links = {}
    for edge in edges_p + edges_p {
        for pair in [{ x: edge.a, y: edge.b }, { x: edge.b, y: edge.a }] {
            if contains(keys(self_links), str(pair.x)) {
                self_links[str(pair.x)] = push(self_links[str(pair.x)], pair.y)
            } else {
                self_links[str(pair.x)] = [pair.y]
            }
        }
    }
    let self_seen = {}
    let self_cycles = 0
    for node in keys(self_links) {
        if contains(keys(self_seen), node) == false {
            self_cycles = self_cycles + 1
            let frontier = [node]
            while len(frontier) > 0 {
                let at = frontier[0]
                frontier = slice(frontier, 1, len(frontier))
                if contains(keys(self_seen), at) == false {
                    self_seen[at] = true
                    for neighbour in self_links[at] {
                        if contains(keys(self_seen), str(neighbour)) == false {
                            frontier = push(frontier, str(neighbour))
                        }
                    }
                }
            }
        }
    }
    assert blocks - self_cycles == 0, "BA6C: a genome is distance 0 from itself"
}

# Rosalind: BA6D — Find a Shortest Transformation of One Genome into Another by 2-Breaks
# https://rosalind.info/problems/ba6d/
#
# Given: Two genomes with circular chromosomes on the same synteny blocks.
# Return: The sequence of genomes along a shortest 2-break transformation.

let genome_p = [[1, -2, -3, 4]]
let genome_q = [[1, 2, -4, -3]]

# BA6C says the distance is blocks minus cycles. This exhibits a path of that
# length, which is the constructive half of the same argument: pick any cycle of
# the combined graph that is not already trivial, and there is always a 2-break
# on P's edges that merges P one step closer to Q while raising the cycle count
# by one. Repeating that reaches Q in exactly blocks-minus-cycles moves.
fn coloured_edges(g) {
    g |> flat_map(|chromosome| {
        let nodes = chromosome |> flat_map(|block|
            if block > 0 then [2 * block - 1, 2 * block] else [0 - 2 * block, 0 - 2 * block - 1])
        range(1, len(chromosome) + 1) |> map(|j| {
            let to_index = if 2 * j == len(nodes) then 0 else 2 * j
            { a: nodes[2 * j - 1], b: nodes[to_index] }
        })
    })
}

fn same_edge(edge, x, y) {
    (edge.a == x and edge.b == y) or (edge.a == y and edge.b == x)
}

fn graph_to_genome(edge_list) {
    let links = {}
    for edge in edge_list {
        links[str(edge.a)] = edge.b
        links[str(edge.b)] = edge.a
    }
    let visited = {}
    let genome = []
    for edge in edge_list {
        if contains(keys(visited), str(edge.a)) == false {
            let cycle = []
            let node = edge.a
            let walking = true
            while walking {
                visited[str(node)] = true
                let partner = links[str(node)]
                visited[str(partner)] = true
                cycle = push(cycle, partner)
                let next_node = if partner % 2 == 1 then partner + 1 else partner - 1
                if contains(keys(visited), str(next_node)) {
                    walking = false
                } else {
                    node = next_node
                }
            }
            let blocks = cycle |> map(|tail|
                if tail % 2 == 1 then (tail + 1) / 2 else 0 - tail / 2)
            let lowest = blocks |> map(|b| abs(b)) |> min()
            let at = (range(0, len(blocks)) |> filter(|i| abs(blocks[i]) == lowest))[0]
            genome = push(genome, range(0, len(blocks)) |> map(|i| blocks[(i + at) % len(blocks)]))
        }
    }
    genome
}

fn written(g) {
    g |> map(|c| "(" + (c |> map(|v| if v > 0 then "+" + str(v) else str(v)) |> join(" ")) + ")")
      |> join("")
}

let red = coloured_edges(genome_p)
let blue = coloured_edges(genome_q)
let steps = [written(graph_to_genome(red))]

# Each round: find a blue edge whose endpoints are joined differently in red, and
# rewire red to match it.
let working = red
let rounds = 0
while rounds < 20 {
    rounds = rounds + 1
    let wrong = blue |> filter(|b| (working |> count_if(|r| same_edge(r, b.a, b.b))) == 0)
    if len(wrong) == 0 { rounds = 100 }
    if rounds < 100 {
        let target = wrong[0]
        # The two red edges currently holding target's endpoints.
        let holding_a = (working |> filter(|r| r.a == target.a or r.b == target.a))[0]
        let holding_b = (working |> filter(|r| r.a == target.b or r.b == target.b))[0]
        let other_a = if holding_a.a == target.a then holding_a.b else holding_a.a
        let other_b = if holding_b.a == target.b then holding_b.b else holding_b.a

        working = (working |> filter(|r|
            same_edge(r, holding_a.a, holding_a.b) == false
            and same_edge(r, holding_b.a, holding_b.b) == false))
            + [{ a: target.a, b: target.b }, { a: other_a, b: other_b }]
        steps = push(steps, written(graph_to_genome(working)))
    }
}

println("Result:")
for step in steps { println("  " + step) }
println("Expected: (+1 -2 -3 +4) / (+1 -2 -3)(+4) / (+1 -2 -4 -3) / (+1 +2 -4 -3)")

fn test_ba6d_two_break_sorting() {
    # The path has to start at P, end at Q, and take exactly the 2-break distance
    # in steps — which is the whole claim being demonstrated.
    assert steps[0] == written(genome_p), "BA6D: must start at P, got " + steps[0]
    let arrived = graph_to_genome(working)
    let magnitudes = sort(arrived |> flat_map(|c| c |> map(|b| abs(b))))
    assert magnitudes == [1, 2, 3, 4], "BA6D: every block must survive"
    # Q reached, compared up to rotation and direction as in BA6K.
    let rotations = |chromosome| range(0, len(chromosome))
        |> map(|shift| join(range(0, len(chromosome))
            |> map(|i| str(chromosome[(i + shift) % len(chromosome)])), " "))
    let flip = |chromosome| reverse(chromosome) |> map(|b| 0 - b)
    let canonical = |chromosome| min(rotations(chromosome) + rotations(flip(chromosome)))
    assert sort(arrived |> map(|c| canonical(c))) == sort(genome_q |> map(|c| canonical(c))),
        "BA6D: the path ends at " + written(arrived) + ", not at Q"
    # Three 2-breaks, so four genomes listed — matching BA6C's formula.
    assert len(steps) == 4, "BA6D: expected 4 genomes, got " + str(len(steps))
}

# Rosalind: BA6E — Find All Shared k-mers of a Pair of Strings
# https://rosalind.info/problems/ba6e/
#
# Given: An integer k and two strings.
# Return: Every pair of positions (x, y) where the strings share a k-mer, taking
# reverse complements into account.

let k = 3
let first_text = "AAACTCATC"
let second_text = "TTTCAAATC"

# The reverse complement counts because a synteny block conserved between two
# genomes may sit on either strand — an inversion flips a segment, and the same
# genes then read backwards. Ignoring that would make every inverted block look
# like a deletion.
#
# Plotting these pairs gives the dot plot from which rearrangements are read off:
# diagonal runs are conserved blocks, and anti-diagonal runs are inverted ones.
let shared = range(0, len(first_text) - k + 1) |> flat_map(|x| {
    let piece = substr(first_text, x, k)
    let flipped = str(reverse_complement(dna(piece)))
    range(0, len(second_text) - k + 1)
        |> filter(|y| {
            let other = substr(second_text, y, k)
            other == piece or other == flipped
        })
        |> map(|y| { x: x, y: y })
})

println("Result:   " + (shared |> map(|p| "(" + str(p.x) + ", " + str(p.y) + ")") |> join(" ")))
println("Expected: (0, 4) (0, 0) (4, 2) (6, 6)   (in any order)")

fn test_ba6e_shared_kmers() {
    let written = sort(shared |> map(|p| "(" + str(p.x) + ", " + str(p.y) + ")"))
    assert written == sort(["(0, 4)", "(0, 0)", "(4, 2)", "(6, 6)"]),
        "BA6E: got " + join(written, " ")
    # Every pair really shares a k-mer, one way round or the other.
    for pair in shared {
        let piece = substr(first_text, pair.x, k)
        let other = substr(second_text, pair.y, k)
        assert other == piece or other == str(reverse_complement(dna(piece))),
            "BA6E: " + piece + " and " + other + " are not a shared k-mer"
    }
    # (0, 4) is the reverse-complement match — AAA against TTT — which a
    # same-strand search would miss entirely.
    assert substr(first_text, 0, k) == "AAA" and substr(second_text, 4, k) == "AAA",
        "BA6E: position 4 of the second string is a direct match"
    assert substr(second_text, 0, k) == "TTT", "BA6E: and position 0 is its reverse complement"
}

# Rosalind: BA6F — Implement ChromosomeToCycle
# https://rosalind.info/problems/ba6f/
#
# Given: A chromosome of n synteny blocks.
# Return: The 2n node numbers ChromosomeToCycle produces.

let chromosome = [1, -2, -3, 4]

# Every block becomes two nodes, its head and its tail, so orientation stops
# being a sign and becomes a direction of travel. That is the representation
#2-breaks are defined on: a rearrangement cuts two edges and reconnects the four
# loose ends, which is easy to say about nodes and awkward to say about signed
# integers.
let cycle_nodes = chromosome |> flat_map(|block|
    if block > 0 then [2 * block - 1, 2 * block] else [0 - 2 * block, 0 - 2 * block - 1])

println("Result:   (" + (cycle_nodes |> map(|n| str(n)) |> join(" ")) + ")")
println("Expected: (1 2 4 3 6 5 7 8)")

fn test_ba6f_chromosome_to_cycle() {
    assert cycle_nodes == [1, 2, 4, 3, 6, 5, 7, 8], "BA6F: got " + str(cycle_nodes)
    # Two nodes per block, and every number from 1 to 2n used exactly once.
    assert len(cycle_nodes) == 2 * len(chromosome), "BA6F: two cycle_nodes per block"
    assert sort(cycle_nodes) == range(1, 2 * len(chromosome) + 1),
        "BA6F: the cycle_nodes must be 1..2n, each once"
    # A positive block reads head-then-tail; a negative one reads the other way,
    # which is the whole encoding.
    assert cycle_nodes[0] < cycle_nodes[1], "BA6F: +1 runs forwards"
    assert cycle_nodes[2] > cycle_nodes[3], "BA6F: -2 runs backwards"
}

# Rosalind: BA6G — Implement CycleToChromosome
# https://rosalind.info/problems/ba6g/
#
# Given: A sequence of 2n node numbers.
# Return: The chromosome whose cycle they are.

let cycle_nodes = [1, 2, 4, 3, 6, 5, 7, 8]

# The inverse of BA6F. Each pair of nodes is one block, and which of the two
# comes first says which way it reads — so the sign is recovered from the order
# rather than stored.
let chromosome = range(0, len(cycle_nodes) / 2) |> map(|j| {
    let head = cycle_nodes[2 * j]
    let tail = cycle_nodes[2 * j + 1]
    if head < tail then tail / 2 else 0 - head / 2
})

fn written(items) {
    "(" + (items |> map(|v| if v > 0 then "+" + str(v) else str(v)) |> join(" ")) + ")"
}

println("Result:   " + written(chromosome))
println("Expected: (+1 -2 -3 +4)")

fn test_ba6g_cycle_to_chromosome() {
    assert written(chromosome) == "(+1 -2 -3 +4)", "BA6G: got " + written(chromosome)
    # It really inverts BA6F: converting back gives the nodes it started from.
    let round_trip = chromosome |> flat_map(|block|
        if block > 0 then [2 * block - 1, 2 * block] else [0 - 2 * block, 0 - 2 * block - 1])
    assert round_trip == cycle_nodes, "BA6G: the round trip does not return the cycle_nodes"
    assert len(chromosome) == len(cycle_nodes) / 2, "BA6G: one block per pair of cycle_nodes"
}

# Rosalind: BA6H — Implement ColoredEdges
# https://rosalind.info/problems/ba6h/
#
# Given: A genome P.
# Return: The coloured edges of its genome graph.

let chromosomes = [[1, -2, -3], [4, 5, -6]]

# Coloured edges are the ones joining *different* blocks — the adjacencies a
# rearrangement can break. The edges inside a block are fixed, so they carry no
# information and are left out; what remains is exactly the structure two genomes
# can be compared on.
fn chromosome_to_cycle(chromosome) {
    chromosome |> flat_map(|block|
        if block > 0 then [2 * block - 1, 2 * block] else [0 - 2 * block, 0 - 2 * block - 1])
}

let edges_list = chromosomes |> flat_map(|chromosome| {
    let nodes = chromosome_to_cycle(chromosome)
    # The last edge wraps round, because a chromosome here is circular.
    range(1, len(chromosome) + 1) |> map(|j| {
        let to_index = if 2 * j == len(nodes) then 0 else 2 * j
        { a: nodes[2 * j - 1], b: nodes[to_index] }
    })
})

println("Result:   " + (edges_list |> map(|e| "(" + str(e.a) + ", " + str(e.b) + ")") |> join(", ")))
println("Expected: (2, 4), (3, 6), (5, 1), (8, 9), (10, 12), (11, 7)")

fn test_ba6h_coloured_edges() {
    let written = edges_list |> map(|e| "(" + str(e.a) + ", " + str(e.b) + ")") |> join(", ")
    assert written == "(2, 4), (3, 6), (5, 1), (8, 9), (10, 12), (11, 7)",
        "BA6H: got " + written
    # One coloured edge per block, since each block has one outgoing adjacency.
    let block_count = chromosomes |> map(|c| len(c)) |> sum()
    assert len(edges_list) == block_count, "BA6H: one coloured edge per block"
    # Every node appears in at most one coloured edge on each side — the graph is
    # a set of disjoint cycles, which is what makes 2-break distance computable.
    let touched = edges_list |> flat_map(|e| [e.a, e.b])
    assert len(unique(touched)) == len(touched), "BA6H: a node cannot have two coloured edges"
}

# Rosalind: BA6I — Implement GraphToGenome
# https://rosalind.info/problems/ba6i/
#
# Given: The coloured edges of a genome graph.
# Return: The genome.

let coloured = [
    { a: 2, b: 4 }, { a: 3, b: 6 }, { a: 5, b: 1 },
    { a: 7, b: 9 }, { a: 10, b: 12 }, { a: 11, b: 8 },
]

# The inverse of BA6H, and the step that makes 2-breaks usable: a 2-break is
# easy to apply to a set of edges and meaningless as an operation on signed
# integers, so genomes are converted to a graph, rearranged, and converted back.
#
# The edges form disjoint cycles; each cycle is one chromosome. Walking a cycle
# means alternating between the coloured edges given here and the implicit
# black edges inside each block, which is why the walk steps by one node between
# coloured edges.
let neighbours = {}
for edge in coloured {
    neighbours[str(edge.a)] = edge.b
    neighbours[str(edge.b)] = edge.a
}

let visited = {}
let chromosomes = []
for edge in coloured {
    if contains(keys(visited), str(edge.a)) == false {
        # Walk the cycle this edge belongs to, collecting node pairs.
        let cycle = []
        let node = edge.a
        let walking = true
        while walking {
            visited[str(node)] = true
            let partner = neighbours[str(node)]
            visited[str(partner)] = true
            cycle = push(cycle, { head: node, tail: partner })
            # The black edge: from `partner` to the other node of its block.
            let next_node = if partner % 2 == 1 then partner + 1 else partner - 1
            if contains(keys(visited), str(next_node)) { walking = false } else { node = next_node }
        }

        # Each collected pair (tail of one block, head of the next) becomes a
        # block once shifted round by one.
        let blocks = cycle |> map(|pair|
            if pair.tail % 2 == 1 then (pair.tail + 1) / 2 else 0 - pair.tail / 2)

        # A circular chromosome has no distinguished starting block, so the walk
        # can begin anywhere and every rotation is the same chromosome. Rotated
        # here to start at the lowest-numbered block, which is what makes the
        # output comparable to the published one.
        let lowest = blocks |> map(|b| abs(b)) |> min()
        let at = (range(0, len(blocks)) |> filter(|i| abs(blocks[i]) == lowest))[0]
        let rotated = range(0, len(blocks)) |> map(|i| blocks[(i + at) % len(blocks)])
        chromosomes = push(chromosomes, rotated)
    }
}

fn written(items) {
    "(" + (items |> map(|v| if v > 0 then "+" + str(v) else str(v)) |> join(" ")) + ")"
}

let shown = chromosomes |> map(|c| written(c)) |> join("")

println("Result:   " + shown)
println("Expected: (+1 -2 -3)(-4 +5 -6)")

fn test_ba6i_graph_to_genome() {
    assert shown == "(+1 -2 -3)(-4 +5 -6)", "BA6I: got " + shown
    assert len(chromosomes) == 2, "BA6I: the edges form two cycles, so two chromosomes"
    # Rotation-invariance is a real property, not a formatting detail: the walk
    # produced (-2 -3 +1) before rotating, which is the same circular chromosome.
    let first_chromosome = chromosomes[0]
    let rotated_once = range(0, len(first_chromosome))
        |> map(|i| first_chromosome[(i + 1) % len(first_chromosome)])
    assert rotated_once != first_chromosome, "BA6I: a rotation is a different listing"
    assert sort(rotated_once) == sort(first_chromosome), "BA6I: but the same chromosome"
    # Every block from 1 to 6 appears exactly once, either way up.
    let magnitudes = sort(chromosomes |> flat_map(|c| c |> map(|b| abs(b))))
    assert magnitudes == range(1, 7), "BA6I: every block must appear once"
}

# Rosalind: BA6J — Implement 2-BreakOnGenomeGraph
# https://rosalind.info/problems/ba6j/
#
# Given: A genome graph and four nodes i, i', j, j'.
# Return: The graph after removing edges (i, i') and (j, j') and adding (i, j)
# and (i', j').

let edge_set = [
    { a: 2, b: 4 }, { a: 3, b: 8 }, { a: 7, b: 5 }, { a: 6, b: 1 },
]
let breakpoint = { i: 1, i_end: 6, j: 3, j_end: 8 }

# A 2-break is the single operation every rearrangement reduces to: cut two
# adjacencies and rejoin the four loose ends the other way. Reversals,
# translocations, fusions and fissions are all this one move — which is exactly
# why BA6C can put a closed form on the distance.
fn same_edge(edge, x, y) {
    (edge.a == x and edge.b == y) or (edge.a == y and edge.b == x)
}

let kept = edge_set |> filter(|edge|
    same_edge(edge, breakpoint.i, breakpoint.i_end) == false and same_edge(edge, breakpoint.j, breakpoint.j_end) == false)

let rejoined = kept + [{ a: breakpoint.j, b: breakpoint.i }, { a: breakpoint.i_end, b: breakpoint.j_end }]

fn written(edges_list) {
    edges_list |> map(|e| "(" + str(e.a) + ", " + str(e.b) + ")") |> join(", ")
}

println("Result:   " + written(rejoined))
println("Expected: (2, 4), (3, 1), (7, 5), (6, 8)   (in any order)")

fn test_ba6j_two_break_on_graph() {
    # The two edges that were cut are gone, the two new ones are present, and
    # everything else is untouched. Order is not meaningful in an edge set.
    assert len(rejoined) == len(edge_set), "BA6J: a 2-break preserves the edge count"
    assert (rejoined |> count_if(|e| same_edge(e, 1, 6))) == 0, "BA6J: (1, 6) should be gone"
    assert (rejoined |> count_if(|e| same_edge(e, 3, 8))) == 0, "BA6J: (3, 8) should be gone"
    assert (rejoined |> count_if(|e| same_edge(e, 3, 1))) == 1, "BA6J: (3, 1) should be present"
    assert (rejoined |> count_if(|e| same_edge(e, 6, 8))) == 1, "BA6J: (6, 8) should be present"
    assert (rejoined |> count_if(|e| same_edge(e, 2, 4))) == 1, "BA6J: (2, 4) is untouched"
    assert (rejoined |> count_if(|e| same_edge(e, 7, 5))) == 1, "BA6J: (7, 5) is untouched"
    # Every node still has exactly one coloured edge, which is what keeps the
    # graph a set of disjoint cycles.
    let touched = rejoined |> flat_map(|e| [e.a, e.b])
    assert len(unique(touched)) == len(touched), "BA6J: a node cannot gain a second edge"
}

# Rosalind: BA6K — Implement 2-BreakOnGenome
# https://rosalind.info/problems/ba6k/
#
# Given: A genome P and four nodes i, i', j, j'.
# Return: The genome after the 2-break.

let chromosomes = [[1, -2, -4, 3]]
let breakpoint = { i: 1, i_end: 6, j: 3, j_end: 8 }

# The three previous problems assembled. A 2-break is meaningless applied to
# signed integers directly, so the genome is converted to a graph (BA6H), the
# break applied there (BA6J), and the result converted back (BA6I). Here it
# splits one chromosome into two, which is a fission — and the same operation
# with different arguments would fuse, invert or translocate.
fn coloured_edges(g) {
    g |> flat_map(|chromosome| {
        let nodes = chromosome |> flat_map(|block|
            if block > 0 then [2 * block - 1, 2 * block] else [0 - 2 * block, 0 - 2 * block - 1])
        range(1, len(chromosome) + 1) |> map(|j| {
            let to_index = if 2 * j == len(nodes) then 0 else 2 * j
            { a: nodes[2 * j - 1], b: nodes[to_index] }
        })
    })
}

fn same_edge(edge, x, y) {
    (edge.a == x and edge.b == y) or (edge.a == y and edge.b == x)
}

let broken = (coloured_edges(chromosomes) |> filter(|edge|
    same_edge(edge, breakpoint.i, breakpoint.i_end) == false and same_edge(edge, breakpoint.j, breakpoint.j_end) == false))
    + [{ a: breakpoint.i, b: breakpoint.j }, { a: breakpoint.i_end, b: breakpoint.j_end }]

# Back to a genome, exactly as in BA6I.
let links = {}
for edge in broken {
    links[str(edge.a)] = edge.b
    links[str(edge.b)] = edge.a
}

let visited = {}
let rebuilt = []
for edge in broken {
    if contains(keys(visited), str(edge.a)) == false {
        let cycle = []
        let node = edge.a
        let walking = true
        while walking {
            visited[str(node)] = true
            let partner = links[str(node)]
            visited[str(partner)] = true
            cycle = push(cycle, partner)
            let next_node = if partner % 2 == 1 then partner + 1 else partner - 1
            if contains(keys(visited), str(next_node)) { walking = false } else { node = next_node }
        }
        let blocks = cycle |> map(|tail|
            if tail % 2 == 1 then (tail + 1) / 2 else 0 - tail / 2)
        # Circular, so rotate to the lowest-numbered block for a stable listing.
        let lowest = blocks |> map(|b| abs(b)) |> min()
        let at = (range(0, len(blocks)) |> filter(|i| abs(blocks[i]) == lowest))[0]
        rebuilt = push(rebuilt, range(0, len(blocks)) |> map(|i| blocks[(i + at) % len(blocks)]))
    }
}

fn written(items) {
    "(" + (items |> map(|v| if v > 0 then "+" + str(v) else str(v)) |> join(" ")) + ")"
}

let shown = rebuilt |> map(|c| written(c)) |> join(" ")

println("Result:   " + shown)
println("Expected: (+2 -1) (-3 +4)   (up to rotation and which chromosome is listed first)")

fn test_ba6k_two_break_on_genome() {
    # One chromosome became two: this 2-break is a fission.
    assert len(rebuilt) == 2, "BA6K: expected a fission into two chromosomes"
    # Every block survives exactly once, either way up — a 2-break rearranges,
    # it never creates or destroys.
    let magnitudes = sort(rebuilt |> flat_map(|c| c |> map(|b| abs(b))))
    assert magnitudes == [1, 2, 3, 4], "BA6K: every block must survive exactly once"
    # A circular chromosome is the same chromosome under rotation *and* under
    # reading it the other way round, which flips every sign as well as the
    # order. This run produced (+1 -2) where the sample shows (+2 -1) — the same
    # chromosome traversed in the opposite direction. Canonicalising over both
    # symmetries is what makes the comparison meaningful.
    let rotations = |chromosome| range(0, len(chromosome))
        |> map(|shift| join(range(0, len(chromosome))
            |> map(|i| str(chromosome[(i + shift) % len(chromosome)])), " "))
    let flip = |chromosome| reverse(chromosome) |> map(|b| 0 - b)
    let canonical = |chromosome| min(rotations(chromosome) + rotations(flip(chromosome)))

    let mine = sort(rebuilt |> map(|c| canonical(c)))
    let published = sort([canonical([2, -1]), canonical([-3, 4])])
    assert mine == published,
        "BA6K: got " + join(mine, " / ") + " against " + join(published, " / ")
    # And the two really are different listings of the same thing.
    assert canonical([1, -2]) == canonical([2, -1]),
        "BA6K: (+1 -2) and (+2 -1) are one circular chromosome"
}