Sequence

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

# Rosalind: DNA — Counting DNA Nucleotides
# https://rosalind.info/problems/dna/
#
# Given: A DNA string s of length at most 1000 nt.
# Return: Four integers counting A, C, G, T.

let s = dna"AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC"

let counts = base_counts(s)
let result = str(counts.A) + " " + str(counts.C) + " " + str(counts.G) + " " + str(counts.T)

println("Result:   " + result)
println("Expected: 20 12 17 21")

fn test_dna_nucleotide_counts() {
    assert result == "20 12 17 21", "DNA: got '" + result + "'"
}

# Rosalind: RNA — Transcribing DNA into RNA
# https://rosalind.info/problems/rna/
#
# Given: A DNA string t.
# Return: The RNA string u, with every T replaced by U.

let t = dna"GATGGAACTTGACTACGTAAATT"
let u = transcribe(t)

println("Result:   " ++ str(u))
println("Expected: GAUGGAACUUGACUACGUAAAUU")

fn test_rna_transcription() {
    assert str(u) == "GAUGGAACUUGACUACGUAAAUU", "RNA: got " ++ str(u)
}

# Rosalind: REVC — Complementing a Strand of DNA
# https://rosalind.info/problems/revc/
#
# Given: A DNA string s of length at most 1000 bp.
# Return: The reverse complement of s.

let s = dna"AAAACCCGGT"
let rc = reverse_complement(s)

println("Result:   " ++ str(rc))
println("Expected: ACCGGGTTTT")

fn test_revc_reverse_complement() {
    assert str(rc) == "ACCGGGTTTT", "REVC: got " ++ str(rc)
}

# Rosalind: HAMM — Counting Point Mutations
# https://rosalind.info/problems/hamm/
#
# Given: Two DNA strings s and t of equal length.
# Return: The Hamming distance between them.

let s = dna"GAGCCTACTAACGGGAT"
let t = dna"CATCGTAATGACGGCCT"

let distance = hamming_distance(s, t)

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

fn test_hamm_point_mutations() {
    assert distance == 7, "HAMM: got " + str(distance)
}

# Rosalind: GC — Computing GC Content
# https://rosalind.info/problems/gc/
#
# Given: At most 10 DNA strings in FASTA format.
# Return: The ID of the string with the highest GC content, and that content.

let records = [
    { id: "Rosalind_6404", seq: dna"CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAATAATTCTGAGG" },
    { id: "Rosalind_5959", seq: dna"CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTTGTCAGCAGACACGC" },
    { id: "Rosalind_0808", seq: dna"CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT" }
]

let scored = records |> map(|r| { id: r.id, gc: gc_content(r.seq) * 100.0 })
let best = scored |> sort_by(|r| r.gc) |> reverse() |> first()

println("Result:   " + best.id + " " + str(best.gc))
println("Expected: Rosalind_0808 60.919540")

fn test_gc_highest_content() {
    assert best.id == "Rosalind_0808", "GC: picked " + best.id
    assert round(best.gc, 6) == 60.91954, "GC: got " + str(best.gc)
}

# Rosalind: TRAN — Transitions and Transversions
# https://rosalind.info/problems/tran/
#
# Given: Two DNA strings of equal length.
# Return: The transition/transversion ratio.

let s1 = "GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGAAGTACGGGCATCAACCCAGTT"
let s2 = "TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGCGGTACGAGTGTTCCTTTGGGT"

# A transition swaps purine for purine (A<->G) or pyrimidine for pyrimidine
# (C<->T); anything else is a transversion.
fn is_purine(base) { base == "A" or base == "G" }

fn is_transition(a, b) { is_purine(a) == is_purine(b) }

let differing = range(0, len(s1)) |> filter(|i| substr(s1, i, 1) != substr(s2, i, 1))
let transitions = differing |> count_if(|i| is_transition(substr(s1, i, 1), substr(s2, i, 1)))
let transversions = len(differing) - transitions

let ratio = float(transitions) / float(transversions)

println("Result:   " + str(ratio))
println("Expected: 1.21428571429")

fn test_tran_ratio() {
    assert round(ratio, 8) == 1.21428571, "TRAN: got " + str(ratio)
}

# Rosalind: CONS — Consensus and Profile
# https://rosalind.info/problems/cons/
#
# Given: A collection of DNA strings of equal length in FASTA format.
# Return: A consensus string, and the profile matrix of base counts per column.

let strings = [
    "ATCCAGCT", "GGGCAACT", "ATGGATCT", "AAGCAACC",
    "TTGGAACT", "ATGCCATT", "ATGGCACT"
]

let bases = ["A", "C", "G", "T"]
# Held as plain strings: str() on a DNA value renders as "DNA(...)", so substr
# would slice the wrapper rather than the sequence.
let width = len(strings[0])

# profile[base][column] — how often each base appears at each position.
let profile = bases |> map(|b| {
    range(0, width) |> map(|i| strings |> count_if(|s| substr(s, i, 1) == b))
})

# The consensus takes whichever base is commonest in each column.
let consensus_string = range(0, width)
    |> map(|i| {
        let column = range(0, 4) |> map(|b| { base: bases[b], n: profile[b][i] })
        let best = column |> sort_by(|e| e.n) |> reverse() |> first()
        best.base
    })
    |> join("")

println("Result:   " + consensus_string)
println("Expected: ATGCAACT")
range(0, 4) |> each(|b| println("  " + bases[b] + ": " + (profile[b] |> map(|n| str(n)) |> join(" "))))

fn test_cons_consensus_and_profile() {
    assert consensus_string == "ATGCAACT", "CONS: got " + consensus_string
    assert profile[0][0] == 5, "CONS: A count at column 1 was " + str(profile[0][0])
    assert profile[3][0] == 1, "CONS: T count at column 1 was " + str(profile[3][0])
}

# Rosalind: REVP — Locating Restriction Sites
# https://rosalind.info/problems/revp/
#
# Given: A DNA string of length at most 1 kbp.
# Return: The position and length of every reverse palindrome of length
# between 4 and 12, as 1-based positions.

let s = "TCAATGCATGCGGGTCTATATGCAT"

# A reverse palindrome reads the same as its own reverse complement, so a
# window qualifies when it equals reverse_complement(window).
fn is_reverse_palindrome(window) {
    str(reverse_complement(dna(window))) == window
}

let found = range(4, 13) |> flat_map(|width| {
    range(0, len(s) - width + 1)
        |> filter(|i| is_reverse_palindrome(substr(s, i, width)))
        |> map(|i| { position: i + 1, length: width })
}) |> sort_by(|hit| hit.position)

println("Result:   " + str(len(found)) + " palindromes")
found |> each(|hit| println("  " + str(hit.position) + " " + str(hit.length)))
println("Expected: 8 — (4,6) (5,4) (6,6) (7,4) (17,4) (18,4) (20,6) (21,4)")

fn test_revp_reverse_palindromes() {
    assert len(found) == 8, "REVP: got " + str(len(found))
    assert found[0].position == 4 and found[0].length == 6, "REVP: first was " + str(found[0])
}

# Rosalind: CORR — Error Correction in Reads
# https://rosalind.info/problems/corr/
#
# Given: A collection of reads of equal length. Each read either appears in the
# dataset at least twice (counting reverse complements), or is a single-symbol
# error away from exactly one such read.
# Return: Every correction, as "wrong->right".

let reads = ["TCATC", "TTCAT", "TCATC", "TGAAA", "GAGGA", "TTTCA", "ATCAA", "TTGAT", "TTTCC"]

fn revcomp(read) { str(reverse_complement(dna(read))) }

fn hamming_of(a, b) {
    range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1))
}

# A read is correct when it, plus its reverse complement, appears at least
# twice across the dataset.
fn occurrences(read) {
    let direct = reads |> count_if(|r| r == read)
    let reversed = reads |> count_if(|r| r == revcomp(read))
    direct + reversed
}

let correct = reads |> filter(|r| occurrences(r) >= 2) |> unique()
let wrong = reads |> filter(|r| occurrences(r) < 2) |> unique()

# Each incorrect read is one substitution from exactly one correct read, in
# either orientation.
let corrections = wrong |> map(|bad| {
    let targets = correct |> flat_map(|good| [good, revcomp(good)])
    let repaired = targets |> filter(|t| hamming_of(bad, t) == 1) |> first()
    bad ++ "->" ++ repaired
})

println("Result:")
corrections |> each(|c| println("  " + c))
println("Expected: TTCAT->TTGAT, GAGGA->GATGA, TTTCC->TTTCA")

fn test_corr_error_correction() {
    assert len(corrections) == 3, "CORR: got " + str(len(corrections))
    assert corrections |> contains("TTCAT->TTGAT"), "CORR: missing TTCAT->TTGAT"
    assert corrections |> contains("GAGGA->GATGA"), "CORR: missing GAGGA->GATGA"
    assert corrections |> contains("TTTCC->TTTCA"), "CORR: missing TTTCC->TTTCA"
}

# Rosalind: KMER — k-Mer Composition
# https://rosalind.info/problems/kmer/
#
# Given: A DNA string s.
# Return: The 4-mer composition of s: the count of every one of the 256
# possible 4-mers, in lexicographic order.

let s = "CTTCGAAAGTTTGGGCCGAGTCTTATAGTCGAATAGCCTCCATATTGGCCTCTATCTGAGCTCTACGCACGCCAGCTCCCGCCCTGGCACTGCTTGATATTGGGTCCTATCTGACCGGACTTGACGTACCTGGGCTGATATTGCCGATCTGATTTGTCCAAAATTGACCATAGTTAAATGGCGGAATGGGCTTTGCCTCGCACCGTAGGCGATTCCCAAGTGCAGGCAAGGATATTGCGAAATATCTTCCAGAGGGCGGAAATTGCTATTCATCCTATTTGGCTTGCTCTTAAGATCACCGATCTAGCTTCGAATCGCCATGCACAGTTTGGGCTCATCGGAAGTCTCGGGCTGGTTTAATTGTGATTGGGCCTTAGCTCTCGGCTGATCTTCGTCACGATATCGCTAGATCGGCTTGACTTCGCATTGGTCTGATTTTTAAGATCTTAGGCTTGACGGCGATCTCGCTCGATCCTGCCTGATAGGCTAGCTCGATTCTCATTAAGCTCTAGCCTTCGGCCTCTGGCTCTTTAAGCGATGCTGCACGATATCTATCGGACTGCTCTTGCTCGATATCGACTGCTTGACGGCTTGATATCGCTCTTGGCTTGCTTGACTCGCTCTTGATCGCTCTAGCTAGCTCGCTCTGATCTTGCACTTGCTCGGCTCTTGCTCGGCTCTTAGCTCGCTCGATCGCTCTAGCTCGCTCGATCTTGCTCGATCTTGCACTGCTCGATCTGCACTTGCTCGATCTTGCTCGATCTAGCT"
let k = 4

let alphabet = ["A", "C", "G", "T"]

# Every 4-mer in lexicographic order, so the counts line up with the ordering
# the problem asks for.
let all_kmers = range(0, k) |> reduce(
    |acc, _| acc |> flat_map(|prefix| alphabet |> map(|b| prefix ++ b)),
    [""]
)

# Count the windows of s once, then look each k-mer up.
let window_list = range(0, len(s) - k + 1) |> map(|i| substr(s, i, k))
let counts = all_kmers |> map(|kmer| window_list |> count_if(|w| w == kmer))

println("Result:   " + str(len(counts)) + " counts, first 12: " + (take(counts, 12) |> map(|c| str(c)) |> join(" ")))
println("Total:    " + str(sum(counts)) + " window_list (expected " + str(len(s) - k + 1) + ")")

fn test_kmer_composition() {
    assert len(counts) == 256, "KMER: got " + str(len(counts)) + " entries"
    assert all_kmers[0] == "AAAA", "KMER: order starts at " + all_kmers[0]
    assert all_kmers[255] == "TTTT", "KMER: order ends at " + all_kmers[255]
    # Every window is counted exactly once, so the counts must total the number
    # of windows — the check that catches an ordering or lookup mistake.
    assert sum(counts) == len(s) - k + 1, "KMER: counts total " + str(sum(counts))
    # Cross-check one entry by counting it independently.
    let aaaa = window_list |> count_if(|w| w == "AAAA")
    assert counts[0] == aaaa, "KMER: AAAA count disagrees"
}