Alignment

NEED — Pairwise Global Alignment

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

Calls NCBI, so it needs a network connection and its answer can change over time.

Reproduces the EMBOSS Needle score of 257 using align() with affine gaps (match 5, mismatch -4, gap_open -9, gap_extend -1). Asserted, but excluded from the hermetic CI gate because it fetches from NCBI.

# Rosalind: NEED — Pairwise Global Alignment
# https://rosalind.info/problems/need/
#
# Given: Two GenBank IDs.
# Return: Maximum global alignment score (DNAfull matrix, gap_open=10, gap_extend=1).

let ids = ["JX205496.1", "JX469991.1"]

# Fetch each record and drop the ">" header line. drop(_, 1) rather than
# tail(): tail() returns the *last* n lines and defaults to 5.
let sequences = ids |> map(|id| drop(split(ncbi_sequence(id), "\n"), 1) |> join(""))

println("Sequences:")
range(0, len(ids)) |> each(|i| println("  " + ids[i] + ": " + str(len(sequences[i])) + " bp"))

# DNAfull scores an ACGT match +5 and a mismatch -4.
#
# EMBOSS charges `gapopen` for the first position of a gap and `gapextend` for
# each further one. This aligner charges `gap_open + gap_extend` to start a gap
# and `gap_extend` to continue it, so gapopen=10 / gapextend=1 becomes
# gap_open = -9 and gap_extend = -1.
let result = align(sequences[0], sequences[1], "global", 5, -4, -1, -9)

println("\nResult:   " + str(result.score))
println("Expected: 257")
println("Match:    " + str(result.score == 257))

println("\nAlignment detail:")
println("  Identity: " + str(result.identity))
println("  Gaps:     " + str(result.gaps))

fn test_need_global_alignment_score() {
    assert result.score == 257, "NEED: expected 257, got " + str(result.score)
}

SUBO — Suboptimal Local Alignment

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

Brute-force O(n^3) window scan; ~25 s. A motif-index approach would be the natural follow-up.

# Rosalind: SUBO — Suboptimal Local Alignment
# https://rosalind.info/problems/subo/
#
# Given: Two DNA strings with a shared inexact repeat (32-40 bp, <=3 changes).
# Return: Total occurrences of the repeat in each string.

let seq1 = "GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTTGACTCTTTTGTTGGCCTTAAATAGATACATATTTGTGCGACTCCACGAGTGATTCGTA"
let seq2 = "ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAACAAGTGTGCACTTAGCCTTGCCGACTCCTTTGTTTGCCTTAAATAGATACATATTTG"

# Strategy: find all 32-bp substrings of seq1 that appear in seq2 with <=3 mismatches
let motif_len = 33

# Helper: count mismatches between two strings of equal length
# (hamming_distance works on sequences; we use manual counting for strings)
let s1 = seq1
let s2 = seq2

# Find the best shared motif by checking each 33-bp window of s1 against all windows of s2
let best_motif = ""
let best_total = 0

let i = 0
while i <= len(s1) - motif_len {
    let candidate = substr(s1, i, motif_len)

    # Count approximate matches in s2
    let hits_s2 = 0
    let j = 0
    while j <= len(s2) - motif_len {
        let target = substr(s2, j, motif_len)
        let dist = hamming_distance(candidate, target)
        if dist <= 3 then
            hits_s2 = hits_s2 + 1
        j = j + 1
    }

    if hits_s2 > 0 then {
        # Also count in s1
        let hits_s1 = 0
        let k = 0
        while k <= len(s1) - motif_len {
            let target = substr(s1, k, motif_len)
            let dist = hamming_distance(candidate, target)
            if dist <= 3 then
                hits_s1 = hits_s1 + 1
            k = k + 1
        }
        let total = hits_s1 + hits_s2
        if total > best_total then {
            best_total = total
            best_motif = candidate
        }
    }
    i = i + 1
}

# Now count non-overlapping occurrences using the best motif
# Count all approximate matches (overlapping) — Rosalind counts overlapping instances
let count1 = 0
let i = 0
while i <= len(s1) - motif_len {
    let window = substr(s1, i, motif_len)
    if hamming_distance(best_motif, window) <= 3 then
        count1 = count1 + 1
    i = i + 1
}

let count2 = 0
let j = 0
while j <= len(s2) - motif_len {
    let window = substr(s2, j, motif_len)
    if hamming_distance(best_motif, window) <= 3 then
        count2 = count2 + 1
    j = j + 1
}

println("Motif:    " + best_motif)
println("Result:   " + str(count1) + " " + str(count2))
println("Expected: 2 2")
println("Match:    " + str(count1 == 2 && count2 == 2))

fn test_subo_repeat_occurrences() {
    assert count1 == 2 && count2 == 2, "SUBO: expected 2 2, got " + str(count1) + " " + str(count2)
}

CLUS — Global Multiple Alignment

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

# Rosalind: CLUS — Global Multiple Alignment
# https://rosalind.info/problems/clus/
#
# Given: Set of DNA strings in FASTA format.
# Return: ID of the string most different from the others.
#
# We compute pairwise edit distances and find the sequence with
# the highest average distance to all others.

let sequences = [
    {id: "Rosalind_18", seq: dna"GACATGTTTGTTTGCCTTAAACTCGTGGCGGCCTAGCCGTAAGTTAAG"},
    {id: "Rosalind_23", seq: dna"ACTCATGTTTGTTTGCCTTAAACTCTTGGCGGCTTAGCCGTAACTTAAG"},
    {id: "Rosalind_51", seq: dna"TCCTATGTTTGTTTGCCTCAAACTCTTGGCGGCCTAGCCGTAAGGTAAG"},
    {id: "Rosalind_7",  seq: dna"CACGTCTGTTCGCCTAAAACTTTGATTGCCGGCCTACGCTAGTTAGTTA"},
    {id: "Rosalind_28", seq: dna"GGGGTCATGGCTGTTTGCCTTAAACCCTTGGCGGCCTAGCCGTAATGTTT"}
]

# Compute pairwise distances and find the most distant sequence
let most_different_id = ""
let max_avg_dist = 0

for i_seq in sequences {
    let total_dist = 0
    let pair_count = 0
    for j_seq in sequences {
        if i_seq.id != j_seq.id then {
            let dist = edit_distance(str(i_seq.seq), str(j_seq.seq))
            total_dist = total_dist + dist
            pair_count = pair_count + 1
        }
    }
    let avg_dist = total_dist / pair_count
    if avg_dist > max_avg_dist then {
        max_avg_dist = avg_dist
        most_different_id = i_seq.id
    }
}

println("Result:   " + most_different_id)
println("Expected: Rosalind_7")
println("Match:    " + str(most_different_id == "Rosalind_7"))

fn test_clus_most_different_sequence() {
    assert most_different_id == "Rosalind_7", "CLUS: expected Rosalind_7, got " + most_different_id
}