Dynamic programming

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

# Rosalind: LGIS — Longest Increasing Subsequence
# https://rosalind.info/problems/lgis/
#
# Given: A positive integer n and a permutation of length n.
# Return: A longest increasing subsequence, then a longest decreasing one.

let permutation = [5, 1, 4, 2, 3]

# O(n^2) dynamic programme: best[i] is the length of the longest run ending at
# i, and prev[i] remembers which index it came from so the run can be rebuilt.
fn longest_run(values, increasing) {
    let n = len(values)
    let best = range(0, n) |> map(|_| 1)
    let prev = range(0, n) |> map(|_| -1)

    let i = 1
    while i < n {
        let j = 0
        while j < i {
            let ordered = if increasing then values[j] < values[i] else values[j] > values[i]
            if ordered and best[j] + 1 > best[i] {
                best = set_at(best, i, best[j] + 1)
                prev = set_at(prev, i, j)
            }
            j = j + 1
        }
        i = i + 1
    }

    let last = 0
    let k = 1
    while k < n {
        if best[k] > best[last] then last = k
        k = k + 1
    }

    let run = []
    let at = last
    while at >= 0 {
        run = concat([values[at]], run)
        at = prev[at]
    }
    run
}

fn set_at(list, index, value) {
    range(0, len(list)) |> map(|i| if i == index then value else list[i])
}

let increasing = longest_run(permutation, true)
let decreasing = longest_run(permutation, false)

println("Increasing: " + (increasing |> map(|x| str(x)) |> join(" ")))
println("Decreasing: " + (decreasing |> map(|x| str(x)) |> join(" ")))
println("Expected:   1 2 3  /  5 4 2")

fn test_lgis_longest_runs() {
    assert len(increasing) == 3, "LGIS: increasing length " + str(len(increasing))
    assert len(decreasing) == 3, "LGIS: decreasing length " + str(len(decreasing))
}

# Rosalind: LCSQ — Finding a Shared Spliced Motif
# https://rosalind.info/problems/lcsq/
#
# Given: Two DNA strings s and t.
# Return: A longest common subsequence of s and t.

let s = "AACCTTGG"
let t = "ACACTGTGA"

# Standard LCS table, then walk it backwards to rebuild one optimal answer.
# Any longest common subsequence is accepted, so the traceback's tie-breaking
# does not matter.
fn lcs_of(a, b) {
    let rows = len(a) + 1
    let cols = len(b) + 1
    let table = range(0, rows) |> map(|_| range(0, cols) |> map(|_| 0))

    let i = 1
    while i < rows {
        let ai = substr(a, i - 1, 1)
        let row = [0]
        let j = 1
        while j < cols {
            let value = if ai == substr(b, j - 1, 1) {
                table[i - 1][j - 1] + 1
            } else {
                max([table[i - 1][j], row[j - 1]])
            }
            row = push(row, value)
            j = j + 1
        }
        table = set_row(table, i, row)
        i = i + 1
    }

    # Traceback from the bottom-right corner.
    let result = ""
    let x = len(a)
    let y = len(b)
    while x > 0 and y > 0 {
        if substr(a, x - 1, 1) == substr(b, y - 1, 1) {
            result = substr(a, x - 1, 1) ++ result
            x = x - 1
            y = y - 1
        } else {
            if table[x - 1][y] >= table[x][y - 1] then x = x - 1 else y = y - 1
        }
    }
    result
}

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

fn is_subsequence(needle, haystack) {
    let at = 0
    let i = 0
    while i < len(haystack) and at < len(needle) {
        if substr(haystack, i, 1) == substr(needle, at, 1) then at = at + 1
        i = i + 1
    }
    at == len(needle)
}

let motif = lcs_of(s, t)

println("Result:   " + motif + " (length " + str(len(motif)) + ")")
println("Expected: a length-6 common subsequence, such as AACTGG")

fn test_lcsq_shared_spliced_motif() {
    assert len(motif) == 6, "LCSQ: length " + str(len(motif)) + " for '" + motif + "'"
    assert is_subsequence(motif, s), "LCSQ: '" + motif + "' is not a subsequence of s"
    assert is_subsequence(motif, t), "LCSQ: '" + motif + "' is not a subsequence of t"
}

# Rosalind: SCSP — Interleaving Two Motifs
# https://rosalind.info/problems/scsp/
#
# Given: Two DNA strings s and t.
# Return: A shortest common supersequence of s and t.

let s = "ATCTGAT"
let t = "TGCATA"

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

# The shortest common supersequence is built from the longest common
# subsequence: walk both strings, emitting shared characters once and the rest
# as they come.
fn lcs_table(a, b) {
    let grid = range(0, len(a) + 1) |> map(|_| range(0, len(b) + 1) |> map(|_| 0))
    let i = 1
    while i <= len(a) {
        let ai = substr(a, i - 1, 1)
        let row = [0]
        let j = 1
        while j <= len(b) {
            let value = if ai == substr(b, j - 1, 1) {
                grid[i - 1][j - 1] + 1
            } else {
                max([grid[i - 1][j], row[j - 1]])
            }
            row = push(row, value)
            j = j + 1
        }
        grid = set_row(grid, i, row)
        i = i + 1
    }
    grid
}

let grid = lcs_table(s, t)

let result = ""
let x = len(s)
let y = len(t)
while x > 0 and y > 0 {
    if substr(s, x - 1, 1) == substr(t, y - 1, 1) {
        result = substr(s, x - 1, 1) ++ result
        x = x - 1
        y = y - 1
    } else {
        if grid[x - 1][y] >= grid[x][y - 1] {
            result = substr(s, x - 1, 1) ++ result
            x = x - 1
        } else {
            result = substr(t, y - 1, 1) ++ result
            y = y - 1
        }
    }
}
result = substr(s, 0, x) ++ substr(t, 0, y) ++ result

fn is_subsequence(needle, haystack) {
    let at = 0
    let i = 0
    while i < len(haystack) and at < len(needle) {
        if substr(haystack, i, 1) == substr(needle, at, 1) then at = at + 1
        i = i + 1
    }
    at == len(needle)
}

println("Result:   " + result + " (length " + str(len(result)) + ")")
println("Expected: a length-9 supersequence, such as ATGCATGAT")

fn test_scsp_shortest_common_supersequence() {
    assert len(result) == 9, "SCSP: length " + str(len(result)) + " for '" + result + "'"
    assert is_subsequence(s, result), "SCSP: s is not a subsequence of the result"
    assert is_subsequence(t, result), "SCSP: t is not a subsequence of the result"
}