Alignment

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

# Rosalind: EDIT — Edit Distance
# https://rosalind.info/problems/edit/
#
# Given: Two protein strings s and t.
# Return: The edit distance between them.

let s = "PLEASANTLY"
let t = "MEANLY"

let distance = edit_distance(s, t)

println("Result:   " + str(distance))
println("Expected: 5")

fn test_edit_distance() {
    assert distance == 5, "EDIT: got " + str(distance)
}

# Rosalind: GLOB — Global Alignment with Scoring Matrix
# https://rosalind.info/problems/glob/
#
# Given: Two protein strings.
# Return: The maximum global alignment score using BLOSUM62 and a constant
# gap penalty of 5.

let s = "PLEASANTLY"
let t = "MEANLY"
let gap = -5.0

let blosum = score_matrix("blosum62")
let residues = blosum.row_names
let width = blosum.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# Needleman-Wunsch with a linear gap penalty. Rows walk s, columns walk t.
fn global_score(a, b, mat, names, cols, gap_penalty) {
    let previous = range(0, len(b) + 1) |> map(|j| float(j) * gap_penalty)
    let i = 1
    while i <= len(a) {
        let current = [float(i) * gap_penalty]
        let ai = substr(a, i - 1, 1)
        let j = 1
        while j <= len(b) {
            let diagonal = previous[j - 1] + substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let up = previous[j] + gap_penalty
            let left = current[j - 1] + gap_penalty
            current = push(current, max([diagonal, up, left]))
            j = j + 1
        }
        previous = current
        i = i + 1
    }
    previous[len(b)]
}

let result = global_score(s, t, blosum, residues, width, gap)

println("Result:   " + str(int(result)))
println("Expected: 8")

fn test_glob_blosum62_global_alignment() {
    assert int(result) == 8, "GLOB: got " + str(result)
}

# Rosalind: LOCA — Local Alignment with Scoring Matrix
# https://rosalind.info/problems/loca/
#
# Given: Two protein strings.
# Return: The maximum local alignment score using PAM250 and a constant gap
# penalty of 5.

let s = "MEANLYPRTEINSTRING"
let t = "PLEASANTLYEINSTEIN"
let gap = -5.0

let pam = score_matrix("pam250")
let residues = pam.row_names
let width = pam.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# Smith-Waterman: identical recurrence to the global case except that a cell
# never drops below zero, and the answer is the best cell anywhere.
fn local_score(a, b, mat, names, cols, gap_penalty) {
    let previous = range(0, len(b) + 1) |> map(|j| 0.0)
    let best = 0.0
    let i = 1
    while i <= len(a) {
        let current = [0.0]
        let ai = substr(a, i - 1, 1)
        let j = 1
        while j <= len(b) {
            let diagonal = previous[j - 1] + substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let cell = max([0.0, diagonal, previous[j] + gap_penalty, current[j - 1] + gap_penalty])
            current = push(current, cell)
            if cell > best then best = cell
            j = j + 1
        }
        previous = current
        i = i + 1
    }
    best
}

let result = local_score(s, t, pam, residues, width, gap)

println("Result:   " + str(int(result)))
println("Expected: 23")

fn test_loca_pam250_local_alignment() {
    assert int(result) == 23, "LOCA: got " + str(result)
}

# Rosalind: GAFF — Global Alignment with Scoring Matrix and Affine Gap Penalty
# https://rosalind.info/problems/gaff/
#
# Given: Two protein strings.
# Return: The maximum global alignment score using BLOSUM62, a gap opening
# penalty of 11 and a gap extension penalty of 1.

let s = "PRTEINS"
let t = "PRTWPSEIN"
let gap_open = -11.0
let gap_extend = -1.0

let blosum = score_matrix("blosum62")
let residues = blosum.row_names
let width = blosum.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# Gotoh's three matrices, held one row at a time:
#   M — s[i] aligned to t[j]
#   X — a gap in t (consuming s)
#   Y — a gap in s (consuming t)
# Opening a gap costs gap_open, each further residue costs gap_extend.
fn affine_score(a, b, mat, names, cols, open_penalty, extend_penalty) {
    let n = len(b)
    let very_low = -1000000.0

    let prev_m = concat([0.0], range(1, n + 1) |> map(|_| very_low))
    let prev_x = concat([very_low], range(1, n + 1) |> map(|_| very_low))
    let prev_y = concat([very_low], range(1, n + 1) |> map(|j| open_penalty + float(j - 1) * extend_penalty))

    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        let row_m = [very_low]
        let row_x = [open_penalty + float(i - 1) * extend_penalty]
        let row_y = [very_low]

        let j = 1
        while j <= n {
            let sub = substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let best_prev = max([prev_m[j - 1], prev_x[j - 1], prev_y[j - 1]])
            row_m = push(row_m, best_prev + sub)
            row_x = push(row_x, max([prev_m[j] + open_penalty, prev_x[j] + extend_penalty]))
            row_y = push(row_y, max([row_m[j - 1] + open_penalty, row_y[j - 1] + extend_penalty]))
            j = j + 1
        }

        prev_m = row_m
        prev_x = row_x
        prev_y = row_y
        i = i + 1
    }

    max([prev_m[n], prev_x[n], prev_y[n]])
}

let result = affine_score(s, t, blosum, residues, width, gap_open, gap_extend)

println("Result:   " + str(int(result)))
println("Expected: 8")

fn test_gaff_affine_global_alignment() {
    assert int(result) == 8, "GAFF: got " + str(result)
}

# Rosalind: OAP — Overlap Alignment
# https://rosalind.info/problems/oap/
#
# Given: Two DNA strings s and t.
# Return: The score of an optimal overlap alignment of s and t, where a suffix
# of s is aligned against a prefix of t. Match +1, mismatch and gap -2.

let s = "CTAAGGGATTCCGGTAATTAGACAG"
let t = "ATAGACCATATGTCAGTGACTGTGTAA"
let match_score = 1.0
let penalty = -2.0

# An overlap alignment is a global alignment with two changes: starting
# anywhere in s is free (first column zero), and ending anywhere in t is free
# (the answer is the best value in the last row).
fn overlap_score(a, b, match_value, mismatch_value) {
    let n = len(b)
    let previous = range(0, n + 1) |> map(|j| float(j) * mismatch_value)

    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        # Free start in s: no penalty accumulated down the first column.
        let current = [0.0]
        let j = 1
        while j <= n {
            let same = ai == substr(b, j - 1, 1)
            let diagonal = previous[j - 1] + (if same then match_value else mismatch_value)
            current = push(current, max([diagonal, previous[j] + mismatch_value, current[j - 1] + mismatch_value]))
            j = j + 1
        }
        previous = current
        i = i + 1
    }

    # Free end in t: the best score anywhere along the final row.
    max(previous)
}

let result = overlap_score(s, t, match_score, penalty)

println("Result:   " + str(int(result)))
println("Expected: 1")

fn test_oap_overlap_alignment() {
    assert int(result) == 1, "OAP: got " + str(result)
}

# Rosalind: SIMS — Finding a Motif with Modifications
# https://rosalind.info/problems/sims/
#
# Given: A DNA string s and a shorter motif t.
# Return: The maximum alignment score of t against any substring of s, with
# match +1 and mismatch/gap -1, plus a substring achieving it.

let s = "GCAAACCATAAGCCCTACGTGCCGCCTGTTTAAACTCGCGAACTGAATCTTCTGCTTCACGGTGAAAGTACCACAATGGTATCACACCCCAAGGAAAC"
let t = "GCCGTCAGGCTGGTGTCCG"

# A fitting alignment: t must be used in full, but the alignment may start and
# end anywhere in s. Free starts come from a zero first row, and the answer is
# the best value in the last row.
let rows = len(t) + 1
let cols = len(s) + 1

# table[i][j] — best score aligning the first i symbols of t ending at j in s.
# Two rows are enough, and each cell is written in place — rebuilding the row
# with push() made this quadratic in list operations.
let previous = range(0, cols) |> map(|_| 0)
let current = range(0, cols) |> map(|_| 0)

let i = 1
while i < rows {
    let ti = substr(t, i - 1, 1)
    current[0] = 0 - i
    let j = 1
    while j < cols {
        let same = ti == substr(s, j - 1, 1)
        let diagonal = previous[j - 1] + (if same then 1 else -1)
        let up = previous[j] - 1
        let left = current[j - 1] - 1
        current[j] = max([diagonal, up, left])
        j = j + 1
    }
    let swap = previous
    previous = current
    current = swap
    i = i + 1
}
let last_value = previous

let result = max(last_value)

println("Result:   " + str(result))
println("Expected: the best fitting-alignment score of t within s")

fn test_sims_motif_with_modifications() {
    # The score cannot exceed a perfect match of t, nor fall below aligning
    # every symbol as a mismatch.
    assert result <= len(t), "SIMS: score " + str(result) + " exceeds |t|"
    assert result >= 0 - len(t), "SIMS: score below the all-mismatch floor"
    # A fitting alignment of t against the substring it selects must reproduce
    # the same score, so recompute it independently at the winning column.
    let best_at = (range(0, len(last_value)) |> filter(|c| last_value[c] == result))[0]
    assert best_at > 0, "SIMS: best column is the empty prefix"
    assert last_value[best_at] == result, "SIMS: winning cell disagrees"
}

# Rosalind: EDTA — Edit Distance Alignment
# https://rosalind.info/problems/edta/
#
# Given: Two protein strings s and t.
# Return: The edit distance between them, together with an optimal alignment.

let s = "PRETTY"
let t = "PRTTEIN"

fn set_row(grid, index, row) {
    range(0, len(grid)) |> map(|i| if i == index then row else grid[i])
}

# Levenshtein table: substitution, deletion and insertion each cost one.
let rows = len(s) + 1
let cols = len(t) + 1
let grid = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 then j else i))

let i = 1
while i < rows {
    let si = substr(s, i - 1, 1)
    let row = [i]
    let j = 1
    while j < cols {
        let cost = if si == substr(t, j - 1, 1) then 0 else 1
        row = push(row, min([grid[i - 1][j - 1] + cost, grid[i - 1][j] + 1, row[j - 1] + 1]))
        j = j + 1
    }
    grid = set_row(grid, i, row)
    i = i + 1
}

let distance = grid[rows - 1][cols - 1]

# Walk the table backwards to recover one alignment. Gaps are written as "-".
let top = ""
let bottom = ""
let x = len(s)
let y = len(t)
while x > 0 or y > 0 {
    let sx = if x > 0 then substr(s, x - 1, 1) else ""
    let ty = if y > 0 then substr(t, y - 1, 1) else ""
    let cost = if x > 0 and y > 0 and sx == ty then 0 else 1
    if x > 0 and y > 0 and grid[x][y] == grid[x - 1][y - 1] + cost {
        top = sx ++ top
        bottom = ty ++ bottom
        x = x - 1
        y = y - 1
    } else {
        if x > 0 and grid[x][y] == grid[x - 1][y] + 1 {
            top = sx ++ top
            bottom = "-" ++ bottom
            x = x - 1
        } else {
            top = "-" ++ top
            bottom = ty ++ bottom
            y = y - 1
        }
    }
}

fn without_gaps(text) {
    range(0, len(text)) |> filter(|k| substr(text, k, 1) != "-") |> map(|k| substr(text, k, 1)) |> join("")
}

let mismatches = range(0, len(top)) |> count_if(|k| substr(top, k, 1) != substr(bottom, k, 1))

println("Result:   " + str(distance))
println("  " + top)
println("  " + bottom)
println("Expected: 4 — PRETTY vs PRTTEIN")

fn test_edta_edit_distance_alignment() {
    assert distance == 4, "EDTA: distance " + str(distance)
    # The alignment must have equal length rows that recover the inputs once
    # gaps are removed, and cost exactly the reported distance.
    assert len(top) == len(bottom), "EDTA: rows differ in length"
    assert without_gaps(top) == s, "EDTA: top row does not spell s"
    assert without_gaps(bottom) == t, "EDTA: bottom row does not spell t"
    assert mismatches == distance, "EDTA: alignment costs " + str(mismatches) + ", not " + str(distance)
}

# Rosalind: CTEA — Counting Optimal Alignments
# https://rosalind.info/problems/ctea/
#
# Given: Two protein strings s and t.
# Return: The number of optimal alignments between them, modulo 134,217,727.

let s = "PLEASANTLY"
let t = "MEANLY"
let modulus = 134217727

fn set_row(table, index, row) {
    range(0, len(table)) |> map(|i| if i == index then row else table[i])
}

let rows = len(s) + 1
let cols = len(t) + 1

# First the edit distances, exactly as in EDTA.
let cost = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 then j else i))
let i = 1
while i < rows {
    let si = substr(s, i - 1, 1)
    let row = [i]
    let j = 1
    while j < cols {
        let step = if si == substr(t, j - 1, 1) then 0 else 1
        row = push(row, min([cost[i - 1][j - 1] + step, cost[i - 1][j] + 1, row[j - 1] + 1]))
        j = j + 1
    }
    cost = set_row(cost, i, row)
    i = i + 1
}

# Then count the paths that achieve those distances. A cell's count is the sum
# of the counts of every predecessor that reaches it at optimal cost, so ties
# multiply out rather than being collapsed.
let ways = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 or j == 0 then 1 else 0))
let a = 1
while a < rows {
    let sa = substr(s, a - 1, 1)
    let row = [1]
    let b = 1
    while b < cols {
        let step = if sa == substr(t, b - 1, 1) then 0 else 1
        let best = cost[a][b]
        # Bound separately: a continued expression cannot start a line with `+`.
        let via_diagonal = if cost[a - 1][b - 1] + step == best then ways[a - 1][b - 1] else 0
        let via_up = if cost[a - 1][b] + 1 == best then ways[a - 1][b] else 0
        let via_left = if cost[a][b - 1] + 1 == best then row[b - 1] else 0
        let total = via_diagonal + via_up + via_left
        row = push(row, total % modulus)
        b = b + 1
    }
    ways = set_row(ways, a, row)
    a = a + 1
}

let distance = cost[rows - 1][cols - 1]
let result = ways[rows - 1][cols - 1]

println("Edit distance: " + str(distance))
println("Result:   " + str(result))
println("Expected: 4 optimal alignments at distance 5")

fn test_ctea_optimal_alignment_count() {
    assert distance == 5, "CTEA: distance " + str(distance)
    assert result == 4, "CTEA: got " + str(result)
    # There is always at least one optimal alignment.
    assert result >= 1, "CTEA: counted none"
}

# Rosalind: GCON — Global Alignment with Constant Gap Penalty
# https://rosalind.info/problems/gcon/
#
# Given: Two protein strings.
# Return: The maximum global alignment score using BLOSUM62 and a constant gap
# penalty of 5 — a run of gaps costs 5 however long it is.

let s = "PLEASANTLY"
let t = "MEANLY"
let gap_open = -5.0
let gap_extend = 0.0

let blosum = score_matrix("blosum62")
let residues = blosum.row_names
let width = blosum.ncol

fn residue_index(names, ch) {
    (range(0, len(names)) |> filter(|i| names[i] == ch))[0]
}

fn substitution(mat, names, cols, a, b) {
    mat.data[residue_index(names, a) * cols + residue_index(names, b)]
}

# The affine recurrence with a zero extension cost: opening a gap is charged
# once and continuing it is free, which is exactly "constant".
fn constant_gap_score(a, b, mat, names, cols, open_penalty, extend_penalty) {
    let n = len(b)
    let very_low = -1000000.0

    let prev_m = concat([0.0], range(1, n + 1) |> map(|_| very_low))
    let prev_x = concat([very_low], range(1, n + 1) |> map(|_| very_low))
    let prev_y = concat([very_low], range(1, n + 1) |> map(|_| open_penalty))

    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        let row_m = [very_low]
        let row_x = [open_penalty]
        let row_y = [very_low]

        let j = 1
        while j <= n {
            let sub = substitution(mat, names, cols, ai, substr(b, j - 1, 1))
            let best_prev = max([prev_m[j - 1], prev_x[j - 1], prev_y[j - 1]])
            row_m = push(row_m, best_prev + sub)
            row_x = push(row_x, max([prev_m[j] + open_penalty, prev_x[j] + extend_penalty]))
            row_y = push(row_y, max([row_m[j - 1] + open_penalty, row_y[j - 1] + extend_penalty]))
            j = j + 1
        }

        prev_m = row_m
        prev_x = row_x
        prev_y = row_y
        i = i + 1
    }

    max([prev_m[n], prev_x[n], prev_y[n]])
}

let result = constant_gap_score(s, t, blosum, residues, width, gap_open, gap_extend)

println("Result:   " + str(int(result)))
println("Expected: 13 — the same pair scores 8 under GLOB's per-residue gap of 5,")
println("          so charging the four-residue gap once instead of four times")
println("          recovers 3 x 5 = 15 minus the 10 already counted.")

fn test_gcon_constant_gap_alignment() {
    assert int(result) == 13, "GCON: got " + str(result)
    # A constant gap penalty can never score worse than the linear one, which
    # charges the same opening plus more for every extra residue.
    assert int(result) >= 8, "GCON: scored below the linear-gap result"
}

# Rosalind: LAFF — Local Alignment with Affine Gap Penalty
# https://rosalind.info/problems/laff/
#
# Given: Two protein strings s and t.
# Return: The maximum local alignment score, and the substrings of s and t that
# achieve it. BLOSUM62, gap opening 11, gap extension 1.

let s = "PLEASANTLY"
let t = "MEANLY"

# Local mode forgives both ends and lets the alignment stop early; the affine
# gap is the gap_open/gap_extend pair. A gap of length L costs 11 + (L-1), so
# opening is -10 on top of the -1 every symbol pays.
let result = align(s, t, "local", 0, 0, -1, -10, "blosum62")

# Local alignment reports the aligned substrings; strip the gaps to recover the
# pieces of s and t themselves.
fn without_gaps(aligned) {
    range(0, len(aligned))
      |> map(|i| substr(aligned, i, 1))
      |> filter(|c| c != "-")
      |> join("")
}

let sub_s = without_gaps(result.aligned_a)
let u = without_gaps(result.aligned_b)

println("Result:   " + str(result.score))
println("          " + sub_s)
println("          " + u)
println("Expected: 12 / LEAS / MEAN")

fn test_laff_local_affine_alignment() {
    assert result.score == 12, "LAFF: got " + str(result.score)
    assert contains(s, sub_s), "LAFF: " + sub_s + " is not a substring of s"
    assert contains(t, u), "LAFF: " + u + " is not a substring of t"
}

# Rosalind: SMGB — Semiglobal Alignment
# https://rosalind.info/problems/smgb/
#
# Given: Two DNA strings s and t.
# Return: The maximum semiglobal alignment score, and an alignment achieving it.
# Match +1, substitution -1, linear gap -1.

let s = dna"CAGCACTTGGATTCTCGG"
let t = dna"CAGCGTGG"

# Semiglobal forgives gaps at both ends of both sequences, so neither overhang
# is charged for — which is what lets a short sequence sit inside a long one
# without paying for the flanks. Global alignment of these two scores far worse.
let result = align(s, t, "semiglobal", 1, -1, -1, 0)
let global_score = align(s, t, "global", 1, -1, -1, 0).score

println("Result:   " + str(result.score))
println("Expected: 4")
println("Alignment:")
println("  " + result.aligned_a)
println("  " + result.aligned_b)
println("(global alignment of the same pair scores " + str(global_score) + ")")

fn test_smgb_semiglobal_alignment() {
    assert result.score == 4, "SMGB: got " + str(result.score)
    # Forgiving the end gaps can only help, never hurt.
    assert result.score >= global_score, "SMGB: semiglobal scored below global"
}

# Rosalind: MGAP — Maximizing the Gap Symbols of an Optimal Alignment
# https://rosalind.info/problems/mgap/
#
# Given: Two DNA strings s and t.
# Return: The largest number of gap symbols that can appear in any maximum-score
# alignment of s and t, for any scoring with match > 0 and both penalties < 0.

let s = dna"AACGTA"
let t = dna"ACACCTA"

# The scoring is left open on purpose, and that is the whole problem: since a
# match is worth something and everything else costs, a maximum-score alignment
# must match as many symbols as it possibly can. That is the longest common
# subsequence. Every symbol outside it is opposite a gap — len(s) - lcs of them
# in one row and len(t) - lcs in the other.
let common = lcs(str(s), str(t))
let gaps = len(str(s)) + len(str(t)) - 2 * len(common)

println("LCS:      " + common)
println("Result:   " + str(gaps))
println("Expected: 3")

fn test_mgap_maximum_gap_symbols() {
    assert gaps == 3, "MGAP: got " + str(gaps)
    # The subsequence has to be a real one of both strings.
    assert is_subsequence(common, str(s)), "MGAP: LCS is not a subsequence of s"
    assert is_subsequence(common, str(t)), "MGAP: LCS is not a subsequence of t"
}

# Rosalind: MULT — Multiple Alignment
# https://rosalind.info/problems/mult/
#
# Given: Four DNA strings of length at most 10.
# Return: A multiple alignment of maximum score, where every mismatched pair in
# a column costs 1 and matched symbols — including two gaps — cost nothing.

let sequences = ["ATATCCG", "TCCG", "ATGTACTG", "ATGTCTG"]

# Two sequences need a table, three a cube, four a hypercube — and the number of
# moves out of each cell is 2^k - 1, every non-empty choice of which sequences
# advance. Fifteen here, against three for a pairwise alignment.
#
# That is the whole lesson of this problem: exact multiple alignment costs
# O(n^k), so it stops being possible at a handful of short sequences. Every tool
# that aligns hundreds of sequences is approximating, usually by aligning pairs
# and merging.
let moves = range(1, 16) |> map(|mask|
    range(0, 4) |> map(|s| (mask / pow(2, s)) % 2 |> floor()))

let lengths = sequences |> map(|s| len(s))

fn key(at) { at |> map(|v| str(v)) |> join(",") }

# A column costs one for every disagreeing pair. Two gaps agree; a gap against a
# base does not.
fn column_cost(symbols) {
    range(0, len(symbols))
        |> flat_map(|i| range(i + 1, len(symbols))
            |> filter(|j| symbols[i] != symbols[j]))
        |> len()
}

let best = {}
let came_from = {}
best[key([0, 0, 0, 0])] = 0

for i in range(0, lengths[0] + 1) {
    for j in range(0, lengths[1] + 1) {
        for k in range(0, lengths[2] + 1) {
            for l in range(0, lengths[3] + 1) {
                let here = [i, j, k, l]
                if i + j + k + l > 0 {
                    let at = key(here)
                    best[at] = 0 - 1000000
                    for move in moves {
                        let previous = range(0, 4) |> map(|s| here[s] - move[s])
                        if (previous |> count_if(|v| v < 0)) == 0 {
                            let symbols = range(0, 4) |> map(|s|
                                if move[s] == 1 then substr(sequences[s], previous[s], 1) else "-")
                            let candidate = best[key(previous)] - column_cost(symbols)
                            if candidate > best[at] {
                                best[at] = candidate
                                came_from[at] = move
                            }
                        }
                    }
                }
            }
        }
    }
}

let score = best[key(lengths)]

let rows = ["", "", "", ""]
let at = lengths
while (at |> sum()) > 0 {
    let move = came_from[key(at)]
    for s in range(0, 4) {
        if move[s] == 1 {
            rows[s] = substr(sequences[s], at[s] - 1, 1) + rows[s]
            at[s] = at[s] - 1
        } else {
            rows[s] = "-" + rows[s]
        }
    }
}

println("Result:   " + str(score))
for line in rows { println("          " + line) }
println("Expected: -18")
println("          ATAT-CCG / -T---CCG / ATGTACTG / ATGT-CTG")

fn test_mult_multiple_alignment() {
    assert score == 0 - 18, "MULT: scored " + str(score)
    # Structural checks, since several alignments reach the optimum.
    let widths = rows |> map(|r| len(r)) |> unique()
    assert len(widths) == 1, "MULT: every row must be the same length"
    for s in range(0, 4) {
        assert replace(rows[s], "-", "") == sequences[s],
            "MULT: row " + str(s) + " is not its original sequence"
    }
    # The alignment shown must actually score what was claimed.
    let recounted = 0 - (range(0, widths[0])
        |> map(|c| column_cost(rows |> map(|r| substr(r, c, 1))))
        |> sum())
    assert recounted == score,
        "MULT: the alignment shown scores " + str(recounted) + ", not " + str(score)
    # No column is entirely gaps — that move does not exist, and one would be
    # free under this scoring.
    let empty_columns = range(0, widths[0])
        |> count_if(|c| (rows |> count_if(|r| substr(r, c, 1) == "-")) == 4)
    assert empty_columns == 0, "MULT: an all-gap column would be free and meaningless"
}

# Rosalind: OSYM — Isolating Symbols in Alignments
# https://rosalind.info/problems/osym/
#
# Given: Two DNA strings.
# Return: The maximum global alignment score, and the sum over all (j,k) of the
# best score of any alignment that pairs s[j] with t[k].
# Matches score +1, mismatches and gaps -1.

let first_string = "ATAGATA"
let second_string = "ACAGGTA"
let gap = 0 - 1

fn substitution(a, b) { if a == b then 1 else 0 - 1 }

# The question is not what the best alignment is, but how good the best
# alignment *through each particular pairing* would be. Recomputing an alignment
# for every one of the 49 pairs would be quadratic work repeated quadratically.
#
# Instead: the best alignment through (j,k) is the best way of reaching it plus
# the best way of leaving it. One table filled forwards and one filled backwards
# answer every pair at once — the same trick BA5K uses to find a middle edge.
fn forward_table(rows, columns, indel) {
    let table = [range(0, len(columns) + 1) |> map(|j| j * indel)]
    for i in range(1, len(rows) + 1) {
        let line = [i * indel]
        for j in range(1, len(columns) + 1) {
            let diagonal = table[i - 1][j - 1]
                + substitution(substr(rows, i - 1, 1), substr(columns, j - 1, 1))
            let up = table[i - 1][j] + indel
            let left = line[j - 1] + indel
            line = push(line, max([diagonal, up, left]))
        }
        table = push(table, line)
    }
    table
}

let forward = forward_table(first_string, second_string, gap)
# The backward table is the forward one on both strings reversed, so only one
# direction has to be written.
let backward_reversed = forward_table(reverse(first_string), reverse(second_string), gap)

let n = len(first_string)
let m = len(second_string)

fn backward_at(table, i, j, rows, columns) {
    # Score of aligning s[i..) with t[j..), read out of the reversed table.
    table[rows - i][columns - j]
}

let best_score = forward[n][m]

let total = range(1, n + 1) |> flat_map(|j| range(1, m + 1) |> map(|k| {
    let paired = substitution(substr(first_string, j - 1, 1), substr(second_string, k - 1, 1))
    forward[j - 1][k - 1] + paired + backward_at(backward_reversed, j, k, n, m)
})) |> sum()

println("Result:   " + str(best_score))
println("          " + str(total))
println("Expected: 3")
println("          -139")

fn test_osym_isolating_symbols() {
    assert best_score == 3, "OSYM: the best alignment scores " + str(best_score)
    assert total == 0 - 139, "OSYM: the sum is " + str(total)

    # No pairing can beat the unconstrained optimum, and at least one must reach
    # it — the best alignment pairs something.
    let every = range(1, n + 1) |> flat_map(|j| range(1, m + 1) |> map(|k| {
        let paired = substitution(substr(first_string, j - 1, 1), substr(second_string, k - 1, 1))
        forward[j - 1][k - 1] + paired + backward_at(backward_reversed, j, k, n, m)
    }))
    assert max(every) == best_score,
        "OSYM: the best constrained score should equal the unconstrained one"
    for value in every {
        assert value <= best_score, "OSYM: a constrained alignment cannot beat the optimum"
    }
    assert len(every) == n * m, "OSYM: one entry per pair of positions"
}