///|
/// alignment_frequencies module example - MSA substitution and column
/// statistics ported from Biopython `Bio.Align`:
/// MultipleSeqAlignment.substitutions (symmetric weighted pair-count
/// matrix) and Alignment.frequencies (per-column letter counts with a
/// '-' gap row), plus observed/expected relative pair frequencies.
///
/// Run: moon run examples/alignment_frequencies_demo
fn main {
println("=== Biopython Bio.Align MSA Statistics Demo ===")
println("Substitution counts and column frequencies of multiple alignments")
println("")
// ---- Section 1: substitution count matrix (Biopython doc example) ----
let rows : Array[String] = ["ACGT", "A--A", "ACGT", "TTTC"]
let sub = @bio.alignment_substitution_matrix(rows) catch {
_ => abort("substitution matrix should build")
}
println("1) Substitution count matrix (MultipleSeqAlignment.substitutions)")
println(" MSA: ACGT / A--A / ACGT / TTTC")
println(sub.to_table())
println(" total pair mass: \{sub.total()}")
println("")
// ---- Section 2: per-column counts and consensus --------------------
let counts = @bio.alignment_letter_counts(rows) catch {
_ => abort("letter counts should build")
}
let freq = counts.column_frequencies()
let gap_k = match counts.letter_index("-") {
Some(i) => i
None => -1
}
println("2) Per-column letter counts (Alignment.frequencies; '-' = gaps)")
println(counts.to_table())
let mut consensus = ""
let mut conserved = 0
let mut j = 0
while j < counts.ncols {
let mut best_k = -1
let mut best_f = 0.0
let mut k = 0
while k < counts.letters.length() {
if k != gap_k && freq[k][j] > best_f {
best_f = freq[k][j]
best_k = k
}
k = k + 1
}
let gap_f = if gap_k >= 0 { freq[gap_k][j] } else { 0.0 }
if best_k >= 0 {
consensus += counts.letters[best_k]
} else {
consensus += "-"
}
if best_f >= 0.75 && gap_f < 0.25 {
conserved = conserved + 1
}
j = j + 1
}
println(" consensus sequence : \{consensus}")
println(
" conserved columns : \{conserved} / \{counts.ncols} (>=75% identity)",
)
println("")
// ---- Section 3: sequence weights ------------------------------------
let wrows : Array[String] = ["GTATC", "AT--C", "CTGTC"]
let wweights : Array[Double] = [0.5, 0.8, 1.0]
let wsub = @bio.alignment_substitution_matrix(wrows, weights=wweights) catch {
_ => abort("weighted substitution matrix should build")
}
println("3) Weighted alignment (sequence weights 0.5 / 0.8 / 1.0)")
println(" MSA: GTATC / AT--C / CTGTC")
println(wsub.to_table())
println(
" total weighted pair mass: \{lod_demo_fmt(wsub.total()).trim()} (unweighted would be 9.000)",
)
println("")
// ---- Section 4: observed vs expected pair frequencies --------------
let bg = sub.background_frequencies()
let rel = sub.relative_frequencies() catch {
_ => abort("relative frequencies should build")
}
println("4) Log-odds scores log2(q_ij / e_ij) with MLE background")
let mut header = " "
let mut li = 0
while li < sub.letters.length() {
header += demo_pad_left(sub.letters[li], 7)
li = li + 1
}
println(header)
let mut a = 0
while a < sub.letters.length() {
let mut line = demo_pad_left(sub.letters[a], 7) + " "
let mut b = 0
while b < sub.letters.length() {
if sub.matrix[a][b] == 0.0 {
line += " ."
} else {
let lod = @math.log2(rel[a][b])
line += lod_demo_fmt(lod)
}
b = b + 1
}
println(line)
a = a + 1
}
let mut bg_line = " background: "
li = 0
while li < sub.letters.length() {
let pv = match bg.get(sub.letters[li]) {
Some(v) => v
None => abort("background should cover all letters")
}
bg_line += sub.letters[li] + "=" + lod_demo_fmt(pv) + " "
li = li + 1
}
println(bg_line)
}
///|
/// Left-pad a string with spaces to the requested width.
fn demo_pad_left(text : String, width : Int32) -> String {
let mut out = text
while out.length() < width {
out = " " + out
}
out
}
///|
/// Fixed-point formatting helper (3 decimals, 7-wide) for the demo.
fn lod_demo_fmt(x : Double) -> String {
let neg = x < 0.0
let v = x.abs() + 0.0005
let whole = v.floor().to_int()
let frac = ((v - v.floor()) * 1000.0).floor().to_int()
let frac_s = if frac >= 100 {
frac.to_string()
} else if frac >= 10 {
"0" + frac.to_string()
} else {
"00" + frac.to_string()
}
let mut s = (if neg { "-" } else { " " }) + whole.to_string() + "." + frac_s
while s.length() < 7 {
s = " " + s
}
s
}