Rna

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

RNAS — Wobble Bonding and RNA Secondary Structures

solvedbrowser + CLI Problem statement Open in the workbench Download .bl

Real RNA pairs U-G nearly as readily as A-U, so a count excluding wobble understates what a molecule can fold into. Memoised with has_key — contains(keys(memo), k) rebuilds the key list on every probe and would make the lookup the bottleneck the memo was meant to remove.

# Rosalind: RNAS — Wobble Bonding and RNA Secondary Structures
# https://rosalind.info/problems/rnas/
#
# Given: An RNA string.
# Return: The number of valid non-crossing matchings of its bonding graph,
# allowing wobble (U-G) pairs, where no bond spans fewer than 4 intervening
# positions.

let rna_string = "AUGCUAGUACGGAGCGAGUCUAGCGAGCGAUGUCGUGAGUACUAUAUAUGCGCAUAAGCCACGU"

# Real RNA does not only pair A-U and C-G. The wobble pair U-G is nearly as
# stable and appears throughout functional RNA, so a structure count that
# excludes it understates what a molecule can fold into.
#
# The other constraint is physical: a hairpin cannot turn in fewer than about
# four bases, so a bond between positions closer than that is impossible however
# well the bases match.
fn can_pair(a, b) {
    (a == "A" and b == "U") or (a == "U" and b == "A")
        or (a == "C" and b == "G") or (a == "G" and b == "C")
        or (a == "U" and b == "G") or (a == "G" and b == "U")
}

# Counting by recursion over intervals: position i either stays unpaired, or
# bonds with some k, which splits the interval into two independent halves. Both
# halves are asked about repeatedly, so without a memo the same interval is
# recomputed exponentially often.
#
# `has_key` is what makes the memo worth having — `contains(keys(memo), k)`
# rebuilds the whole key list on every probe, turning the lookup itself into the
# bottleneck it was meant to remove.
let memo = {}

fn count_matchings(i, j) {
    if j - i < 5 { return 1 }
    let key = str(i) + "," + str(j)
    if has_key(memo, key) { return memo[key] }

    # i unpaired.
    let total = count_matchings(i + 1, j)
    # i bonded to k, which must leave four bases inside the loop.
    for k in range(i + 4, j) {
        if can_pair(substr(rna_string, i, 1), substr(rna_string, k, 1)) {
            total = total + count_matchings(i + 1, k) * count_matchings(k + 1, j)
        }
    }
    memo[key] = total
    total
}

let answer = count_matchings(0, len(rna_string))

println("Result:   " + str(answer))
println("Expected: 284850219977421")

fn test_rnas_wobble_bonding() {
    assert answer == 284850219977421, "RNAS: got " + str(answer)
    # The memo is doing real work: far fewer intervals than the recursion visits.
    assert len(keys(memo)) < len(rna_string) * len(rna_string),
        "RNAS: the memo should hold at most one entry per interval"
    assert len(keys(memo)) > 0, "RNAS: the memo should be used at all"

    # Wobble pairs genuinely matter — U-G is accepted where a strict
    # Watson-Crick rule would reject it.
    assert can_pair("U", "G") and can_pair("G", "U"), "RNAS: wobble pairs bond"
    assert can_pair("A", "G") == false, "RNAS: A-G does not"
    # A string too short to turn a hairpin has exactly one matching: the empty one.
    assert count_matchings(0, 4) == 1, "RNAS: nothing can bond within four bases"
    # The empty matching is always counted, so the total is never zero.
    assert answer > 0, "RNAS: the empty matching always counts"
}