--- title: Rosalind — Bioinformatics Stronghold version: 0.1.0 abstract: Worked solutions to all 105 problems of the Rosalind Bioinformatics Stronghold. --- # Rosalind — Bioinformatics Stronghold Worked solutions to all 105 problems of the Rosalind Bioinformatics Stronghold. Generated from `packs/rosalind-stronghold/pack.toml`. Run the whole notebook with `bl notebook rosalind-stronghold.bln`. ## DNA — Counting DNA Nucleotides [Problem statement](https://rosalind.info/problems/dna/) ```biolang # 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 + "'" } ``` ## RNA — Transcribing DNA into RNA [Problem statement](https://rosalind.info/problems/rna/) ```biolang # 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) } ``` ## REVC — Complementing a Strand of DNA [Problem statement](https://rosalind.info/problems/revc/) ```biolang # 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) } ``` ## HAMM — Counting Point Mutations [Problem statement](https://rosalind.info/problems/hamm/) ```biolang # 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) } ``` ## FIB — Rabbits and Recurrence Relations [Problem statement](https://rosalind.info/problems/fib/) ```biolang # Rosalind: FIB — Rabbits and Recurrence Relations # https://rosalind.info/problems/fib/ # # Given: Positive integers n <= 40 and k <= 5. # Return: The total rabbit pairs after n months, when every pair of # reproduction-age rabbits produces a litter of k pairs each month. let n = 5 let k = 3 fn rabbit_pairs(months, litter) { let previous = 1 let current = 1 let month = 3 while month <= months { let next = current + previous * litter previous = current current = next month = month + 1 } if months <= 2 then 1 else current } let result = rabbit_pairs(n, k) println("Result: " + str(result)) println("Expected: 19") fn test_fib_rabbit_pairs() { assert result == 19, "FIB: got " + str(result) } ``` ## GC — Computing GC Content [Problem statement](https://rosalind.info/problems/gc/) ```biolang # 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) } ``` ## PROT — Translating RNA into Protein [Problem statement](https://rosalind.info/problems/prot/) ```biolang # Rosalind: PROT — Translating RNA into Protein # https://rosalind.info/problems/prot/ # # Given: An RNA string s corresponding to a strand of mRNA. # Return: The protein string encoded by s. let s = rna"AUGGCCAUGGCGCCCAGAACUGAGAUCAAUAGUACCCGUAUUAACGGGUGA" # translate() stops at the first stop codon, which is what the problem asks for. let peptide = translate(s) println("Result: " ++ str(peptide)) println("Expected: MAMAPRTEINSTRING") fn test_prot_translation() { assert str(peptide) == "MAMAPRTEINSTRING", "PROT: got " ++ str(peptide) } ``` ## SUBS — Finding a Motif in DNA [Problem statement](https://rosalind.info/problems/subs/) ```biolang # Rosalind: SUBS — Finding a Motif in DNA # https://rosalind.info/problems/subs/ # # Given: Two DNA strings s and t. # Return: All 1-based locations of t as a substring of s, including overlaps. let s = dna"GATATATGCATATACTT" let t = dna"ATAT" # find_motif reports 0-based positions; Rosalind counts from 1. let positions = find_motif(s, t) |> map(|p| p + 1) let result = positions |> map(|p| str(p)) |> join(" ") println("Result: " + result) println("Expected: 2 4 10") fn test_subs_motif_positions() { assert result == "2 4 10", "SUBS: got '" + result + "'" } ``` ## FIBD — Mortal Fibonacci Rabbits [Problem statement](https://rosalind.info/problems/fibd/) ```biolang # Rosalind: FIBD — Mortal Fibonacci Rabbits # https://rosalind.info/problems/fibd/ # # Given: Positive integers n <= 100 and m <= 20. # Return: The number of rabbit pairs alive after n months, when every pair # lives for exactly m months. let n = 6 let m = 3 # ages[i] = pairs currently i months old. Each month the adults breed and the # oldest cohort dies. fn mortal_pairs(months, lifespan) { let ages = range(0, lifespan) |> map(|i| if i == 0 then 1 else 0) let month = 1 while month < months { let newborns = sum(drop(ages, 1)) ages = concat([newborns], take(ages, lifespan - 1)) month = month + 1 } sum(ages) } let result = mortal_pairs(n, m) println("Result: " + str(result)) println("Expected: 4") fn test_fibd_mortal_rabbits() { assert result == 4, "FIBD: got " + str(result) } ``` ## IPRB — Mendel's First Law [Problem statement](https://rosalind.info/problems/iprb/) ```biolang # Rosalind: IPRB — Mendel's First Law # https://rosalind.info/problems/iprb/ # # Given: Positive integers k (homozygous dominant), m (heterozygous) and # n (homozygous recessive) organisms. # Return: The probability that two randomly chosen mating organisms produce an # individual with a dominant allele. let k = 2.0 let m = 2.0 let n = 2.0 let total = k + m + n # Work out the recessive probability and subtract: only mm, mn and nn pairings # can produce a recessive offspring. # A continued expression keeps its operator at the end of the line. let mm = (m / total) * ((m - 1.0) / (total - 1.0)) * 0.25 let mn = (m / total) * (n / (total - 1.0)) * 0.5 let nm = (n / total) * (m / (total - 1.0)) * 0.5 let nn = (n / total) * ((n - 1.0) / (total - 1.0)) let p_recessive = mm + mn + nm + nn let result = 1.0 - p_recessive println("Result: " + str(round(result, 5))) println("Expected: 0.78333") fn test_iprb_dominant_probability() { assert round(result, 5) == 0.78333, "IPRB: got " + str(result) } ``` ## IEV — Calculating Expected Offspring [Problem statement](https://rosalind.info/problems/iev/) ```biolang # Rosalind: IEV — Calculating Expected Offspring # https://rosalind.info/problems/iev/ # # Given: Six integers, the number of couples with each genotype pairing. # Return: The expected number of offspring displaying the dominant phenotype, # assuming every couple has exactly two offspring. let couples = [1, 0, 0, 1, 0, 1] # Probability a single offspring shows the dominant phenotype, per pairing: # AA-AA, AA-Aa, AA-aa are certain; Aa-Aa is 3/4; Aa-aa is 1/2; aa-aa is 0. let dominant_probability = [1.0, 1.0, 1.0, 0.75, 0.5, 0.0] let offspring_per_couple = 2.0 let result = range(0, 6) |> map(|i| float(couples[i]) * dominant_probability[i] * offspring_per_couple) |> sum() println("Result: " + str(result)) println("Expected: 3.5") fn test_iev_expected_offspring() { assert result == 3.5, "IEV: got " + str(result) } ``` ## MRNA — Inferring mRNA from Protein [Problem statement](https://rosalind.info/problems/mrna/) ```biolang # Rosalind: MRNA — Inferring mRNA from Protein # https://rosalind.info/problems/mrna/ # # Given: A protein string of length at most 1000 aa. # Return: The number of RNA strings that could have produced it, modulo # 1,000,000 — remembering the stop codon. let protein_string = "MA" # Codons per amino acid in the standard genetic code. let codon_counts = { A: 4, C: 2, D: 2, E: 2, F: 2, G: 4, H: 2, I: 3, K: 2, L: 6, M: 1, N: 2, P: 4, Q: 2, R: 6, S: 6, T: 4, V: 4, W: 1, Y: 2 } let stop_codons = 3 let modulus = 1000000 let result = range(0, len(protein_string)) |> map(|i| codon_counts[substr(protein_string, i, 1)]) |> reduce(|acc, n| (acc * n) % modulus, stop_codons) println("Result: " + str(result)) println("Expected: 12") fn test_mrna_possible_rna_strings() { assert result == 12, "MRNA: got " + str(result) } ``` ## PRTM — Calculating Protein Mass [Problem statement](https://rosalind.info/problems/prtm/) ```biolang # Rosalind: PRTM — Calculating Protein Mass # https://rosalind.info/problems/prtm/ # # Given: A protein string P of length at most 1000 aa. # Return: The total weight of P, using the monoisotopic mass table. let protein_string = "SKADYEK" let monoisotopic = { A: 71.03711, C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841, G: 57.02146, H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406, M: 131.04049, N: 114.04293, P: 97.05276, Q: 128.05858, R: 156.10111, S: 87.03203, T: 101.04768, V: 99.06841, W: 186.07931, Y: 163.06333 } let result = range(0, len(protein_string)) |> map(|i| monoisotopic[substr(protein_string, i, 1)]) |> sum() println("Result: " + str(round(result, 3))) println("Expected: 821.392") fn test_prtm_protein_mass() { assert round(result, 3) == 821.392, "PRTM: got " + str(result) } ``` ## PERM — Enumerating Gene Orders [Problem statement](https://rosalind.info/problems/perm/) ```biolang # Rosalind: PERM — Enumerating Gene Orders # https://rosalind.info/problems/perm/ # # Given: A positive integer n <= 7. # Return: The total number of permutations of length n, followed by a list of # all such permutations. let n = 3 fn permutations_of(items) { if len(items) <= 1 { [items] } else { range(0, len(items)) |> flat_map(|i| { let chosen = items[i] let rest = concat(take(items, i), drop(items, i + 1)) permutations_of(rest) |> map(|p| concat([chosen], p)) }) } } let orders = permutations_of(range(1, n + 1)) println("Result: " + str(len(orders))) println("Expected: 6") orders |> each(|p| println(" " + (p |> map(|x| str(x)) |> join(" ")))) fn test_perm_permutation_count() { assert len(orders) == 6, "PERM: got " + str(len(orders)) assert len(unique(orders |> map(|p| join(map(p, |x| str(x)), " ")))) == 6, "PERM: duplicates" } ``` ## PPER — Partial Permutations [Problem statement](https://rosalind.info/problems/pper/) ```biolang # Rosalind: PPER — Partial Permutations # https://rosalind.info/problems/pper/ # # Given: Positive integers n and k with 100 >= n > 0 and 10 >= k > 0. # Return: The number of partial permutations P(n, k), modulo 1,000,000. let n = 21 let k = 7 let modulus = 1000000 # P(n, k) = n * (n-1) * ... * (n-k+1), reduced at every step to stay small. let result = range(0, k) |> reduce(|acc, i| (acc * (n - i)) % modulus, 1) println("Result: " + str(result)) println("Expected: 51200") fn test_pper_partial_permutations() { assert result == 51200, "PPER: got " + str(result) } ``` ## EDIT — Edit Distance [Problem statement](https://rosalind.info/problems/edit/) ```biolang # Rosalind: EDIT — Edit Distance # https://rosalind.info/problems/edit/ # # Given: Two protein strings s and t. # Return: The edit distance between them. let s = "PLEASANTLY" let t = "MEANLY" let distance = edit_distance(s, t) println("Result: " + str(distance)) println("Expected: 5") fn test_edit_distance() { assert distance == 5, "EDIT: got " + str(distance) } ``` ## GLOB — Global Alignment with Scoring Matrix [Problem statement](https://rosalind.info/problems/glob/) ```biolang # Rosalind: GLOB — Global Alignment with Scoring Matrix # https://rosalind.info/problems/glob/ # # Given: Two protein strings. # Return: The maximum global alignment score using BLOSUM62 and a constant # gap penalty of 5. let s = "PLEASANTLY" let t = "MEANLY" let gap = -5.0 let blosum = score_matrix("blosum62") let residues = blosum.row_names let width = blosum.ncol fn residue_index(names, ch) { (range(0, len(names)) |> filter(|i| names[i] == ch))[0] } fn substitution(mat, names, cols, a, b) { mat.data[residue_index(names, a) * cols + residue_index(names, b)] } # Needleman-Wunsch with a linear gap penalty. Rows walk s, columns walk t. fn global_score(a, b, mat, names, cols, gap_penalty) { let previous = range(0, len(b) + 1) |> map(|j| float(j) * gap_penalty) let i = 1 while i <= len(a) { let current = [float(i) * gap_penalty] let ai = substr(a, i - 1, 1) let j = 1 while j <= len(b) { let diagonal = previous[j - 1] + substitution(mat, names, cols, ai, substr(b, j - 1, 1)) let up = previous[j] + gap_penalty let left = current[j - 1] + gap_penalty current = push(current, max([diagonal, up, left])) j = j + 1 } previous = current i = i + 1 } previous[len(b)] } let result = global_score(s, t, blosum, residues, width, gap) println("Result: " + str(int(result))) println("Expected: 8") fn test_glob_blosum62_global_alignment() { assert int(result) == 8, "GLOB: got " + str(result) } ``` ## LOCA — Local Alignment with Scoring Matrix [Problem statement](https://rosalind.info/problems/loca/) ```biolang # Rosalind: LOCA — Local Alignment with Scoring Matrix # https://rosalind.info/problems/loca/ # # Given: Two protein strings. # Return: The maximum local alignment score using PAM250 and a constant gap # penalty of 5. let s = "MEANLYPRTEINSTRING" let t = "PLEASANTLYEINSTEIN" let gap = -5.0 let pam = score_matrix("pam250") let residues = pam.row_names let width = pam.ncol fn residue_index(names, ch) { (range(0, len(names)) |> filter(|i| names[i] == ch))[0] } fn substitution(mat, names, cols, a, b) { mat.data[residue_index(names, a) * cols + residue_index(names, b)] } # Smith-Waterman: identical recurrence to the global case except that a cell # never drops below zero, and the answer is the best cell anywhere. fn local_score(a, b, mat, names, cols, gap_penalty) { let previous = range(0, len(b) + 1) |> map(|j| 0.0) let best = 0.0 let i = 1 while i <= len(a) { let current = [0.0] let ai = substr(a, i - 1, 1) let j = 1 while j <= len(b) { let diagonal = previous[j - 1] + substitution(mat, names, cols, ai, substr(b, j - 1, 1)) let cell = max([0.0, diagonal, previous[j] + gap_penalty, current[j - 1] + gap_penalty]) current = push(current, cell) if cell > best then best = cell j = j + 1 } previous = current i = i + 1 } best } let result = local_score(s, t, pam, residues, width, gap) println("Result: " + str(int(result))) println("Expected: 23") fn test_loca_pam250_local_alignment() { assert int(result) == 23, "LOCA: got " + str(result) } ``` ## GAFF — Global Alignment with Scoring Matrix and Affine Gap Penalty [Problem statement](https://rosalind.info/problems/gaff/) ```biolang # Rosalind: GAFF — Global Alignment with Scoring Matrix and Affine Gap Penalty # https://rosalind.info/problems/gaff/ # # Given: Two protein strings. # Return: The maximum global alignment score using BLOSUM62, a gap opening # penalty of 11 and a gap extension penalty of 1. let s = "PRTEINS" let t = "PRTWPSEIN" let gap_open = -11.0 let gap_extend = -1.0 let blosum = score_matrix("blosum62") let residues = blosum.row_names let width = blosum.ncol fn residue_index(names, ch) { (range(0, len(names)) |> filter(|i| names[i] == ch))[0] } fn substitution(mat, names, cols, a, b) { mat.data[residue_index(names, a) * cols + residue_index(names, b)] } # Gotoh's three matrices, held one row at a time: # M — s[i] aligned to t[j] # X — a gap in t (consuming s) # Y — a gap in s (consuming t) # Opening a gap costs gap_open, each further residue costs gap_extend. fn affine_score(a, b, mat, names, cols, open_penalty, extend_penalty) { let n = len(b) let very_low = -1000000.0 let prev_m = concat([0.0], range(1, n + 1) |> map(|_| very_low)) let prev_x = concat([very_low], range(1, n + 1) |> map(|_| very_low)) let prev_y = concat([very_low], range(1, n + 1) |> map(|j| open_penalty + float(j - 1) * extend_penalty)) let i = 1 while i <= len(a) { let ai = substr(a, i - 1, 1) let row_m = [very_low] let row_x = [open_penalty + float(i - 1) * extend_penalty] let row_y = [very_low] let j = 1 while j <= n { let sub = substitution(mat, names, cols, ai, substr(b, j - 1, 1)) let best_prev = max([prev_m[j - 1], prev_x[j - 1], prev_y[j - 1]]) row_m = push(row_m, best_prev + sub) row_x = push(row_x, max([prev_m[j] + open_penalty, prev_x[j] + extend_penalty])) row_y = push(row_y, max([row_m[j - 1] + open_penalty, row_y[j - 1] + extend_penalty])) j = j + 1 } prev_m = row_m prev_x = row_x prev_y = row_y i = i + 1 } max([prev_m[n], prev_x[n], prev_y[n]]) } let result = affine_score(s, t, blosum, residues, width, gap_open, gap_extend) println("Result: " + str(int(result))) println("Expected: 8") fn test_gaff_affine_global_alignment() { assert int(result) == 8, "GAFF: got " + str(result) } ``` ## OAP — Overlap Alignment [Problem statement](https://rosalind.info/problems/oap/) ```biolang # Rosalind: OAP — Overlap Alignment # https://rosalind.info/problems/oap/ # # Given: Two DNA strings s and t. # Return: The score of an optimal overlap alignment of s and t, where a suffix # of s is aligned against a prefix of t. Match +1, mismatch and gap -2. let s = "CTAAGGGATTCCGGTAATTAGACAG" let t = "ATAGACCATATGTCAGTGACTGTGTAA" let match_score = 1.0 let penalty = -2.0 # An overlap alignment is a global alignment with two changes: starting # anywhere in s is free (first column zero), and ending anywhere in t is free # (the answer is the best value in the last row). fn overlap_score(a, b, match_value, mismatch_value) { let n = len(b) let previous = range(0, n + 1) |> map(|j| float(j) * mismatch_value) let i = 1 while i <= len(a) { let ai = substr(a, i - 1, 1) # Free start in s: no penalty accumulated down the first column. let current = [0.0] let j = 1 while j <= n { let same = ai == substr(b, j - 1, 1) let diagonal = previous[j - 1] + (if same then match_value else mismatch_value) current = push(current, max([diagonal, previous[j] + mismatch_value, current[j - 1] + mismatch_value])) j = j + 1 } previous = current i = i + 1 } # Free end in t: the best score anywhere along the final row. max(previous) } let result = overlap_score(s, t, match_score, penalty) println("Result: " + str(int(result))) println("Expected: 1") fn test_oap_overlap_alignment() { assert int(result) == 1, "OAP: got " + str(result) } ``` ## SSET — Counting Subsets [Problem statement](https://rosalind.info/problems/sset/) ```biolang # Rosalind: SSET — Counting Subsets # https://rosalind.info/problems/sset/ # # Given: A positive integer n <= 1000. # Return: The total number of subsets of {1, 2, ..., n}, modulo 1,000,000. let n = 3 let modulus = 1000000 # 2^n, reduced at every step so the value never grows beyond the modulus. let result = range(0, n) |> reduce(|acc, _| (acc * 2) % modulus, 1) println("Result: " + str(result)) println("Expected: 8") fn test_sset_subset_count() { assert result == 8, "SSET: got " + str(result) } ``` ## PMCH — Perfect Matchings and RNA Secondary Structures [Problem statement](https://rosalind.info/problems/pmch/) ```biolang # Rosalind: PMCH — Perfect Matchings and RNA Secondary Structures # https://rosalind.info/problems/pmch/ # # Given: An RNA string s with the same number of A as U and of G as C. # Return: The total number of perfect matchings of basepair edges. let s = "AGCUAGUCAU" fn count_base(seq, base) { range(0, len(seq)) |> count_if(|i| substr(seq, i, 1) == base) } fn factorial(n) { range(1, n + 1) |> reduce(|acc, i| acc * i, 1) } # Every A pairs with some U and every G with some C, independently, so the # count is |A|! * |G|!. let adenine = count_base(s, "A") let guanine = count_base(s, "G") let result = factorial(adenine) * factorial(guanine) println("Result: " + str(result)) println("Expected: 12") fn test_pmch_perfect_matchings() { assert result == 12, "PMCH: got " + str(result) } ``` ## MMCH — Maximum Matchings and RNA Secondary Structures [Problem statement](https://rosalind.info/problems/mmch/) ```biolang # Rosalind: MMCH — Maximum Matchings and RNA Secondary Structures # https://rosalind.info/problems/mmch/ # # Given: An RNA string s. # Return: The total number of maximum matchings of basepair edges. let s = "AUGCUUC" fn count_base(seq, base) { range(0, len(seq)) |> count_if(|i| substr(seq, i, 1) == base) } # P(n, k) — the number of ways to pair k of the scarcer base with n of the # commoner one. fn partial_permutations(n, k) { range(0, k) |> reduce(|acc, i| acc * (n - i), 1) } fn pairing_ways(first, second) { partial_permutations(max([first, second]), min([first, second])) } # A continued expression keeps its operator at the end of the line, so these are # bound separately rather than split across a leading `*`. let au_ways = pairing_ways(count_base(s, "A"), count_base(s, "U")) let gc_ways = pairing_ways(count_base(s, "G"), count_base(s, "C")) let result = au_ways * gc_ways println("Result: " + str(result)) println("Expected: 6") fn test_mmch_maximum_matchings() { assert result == 6, "MMCH: got " + str(result) } ``` ## CAT — Catalan Numbers and RNA Secondary Structures [Problem statement](https://rosalind.info/problems/cat/) ```biolang # Rosalind: CAT — Catalan Numbers and RNA Secondary Structures # https://rosalind.info/problems/cat/ # # Given: An RNA string s. # Return: The total number of noncrossing perfect matchings of basepair edges, # modulo 1,000,000. let s = "AUAU" let modulus = 1000000 fn complements(a, b) { (a == "A" and b == "U") or (a == "U" and b == "A") or (a == "G" and b == "C") or (a == "C" and b == "G") } # Noncrossing means the first base pairs with some k, splitting the rest into # an inside and an outside that cannot interleave — so the counts multiply. # Only even-length intervals can be matched at all. fn noncrossing(seq, lo, hi) { let width = hi - lo if width <= 0 { 1 } else { let first = substr(seq, lo, 1) let total = range(lo + 1, hi) |> filter(|k| ((k - lo) % 2 == 1) and complements(first, substr(seq, k, 1))) |> map(|k| (noncrossing(seq, lo + 1, k) * noncrossing(seq, k + 1, hi)) % modulus) |> sum() total % modulus } } let result = noncrossing(s, 0, len(s)) println("Result: " + str(result)) println("Expected: 2") fn test_cat_noncrossing_matchings() { assert result == 2, "CAT: got " + str(result) } ``` ## LIA — Independent Alleles [Problem statement](https://rosalind.info/problems/lia/) ```biolang # Rosalind: LIA — Independent Alleles # https://rosalind.info/problems/lia/ # # Given: Integers k <= 7 and N <= 2^k. Tom is AaBb; every organism mates with # an AaBb partner and has two children. # Return: The probability that at least N organisms in generation k are AaBb. let k = 2 let n = 1 # Mendel's second law makes the two loci independent, so each child is AaBb # with probability 1/4 regardless of generation. Generation k holds 2^k # organisms, so the count is Binomial(2^k, 1/4). let population = int(pow(2.0, float(k))) let p = 0.25 fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) } fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) } # P(at least N) = 1 - P(fewer than N). let below = range(0, n) |> map(|i| float(choose(population, i)) * pow(p, float(i)) * pow(1.0 - p, float(population - i))) |> sum() let result = 1.0 - below println("Result: " + str(round(result, 3))) println("Expected: 0.684") fn test_lia_independent_alleles() { assert round(result, 3) == 0.684, "LIA: got " + str(result) } ``` ## SEXL — Sex-Linked Inheritance [Problem statement](https://rosalind.info/problems/sexl/) ```biolang # Rosalind: SEXL — Sex-Linked Inheritance # https://rosalind.info/problems/sexl/ # # Given: An array A of allele frequencies for a gene on the X chromosome. # Return: For each frequency, the probability that a randomly selected female # is a carrier of the recessive allele. let freqs = [0.1, 0.5, 0.8] # Males carry one X, so the recessive allele frequency is the given value. A # female is heterozygous with probability 2q(1-q). let results = freqs |> map(|q| 2.0 * q * (1.0 - q)) let formatted = results |> map(|r| str(round(r, 3))) |> join(" ") println("Result: " + formatted) println("Expected: 0.18 0.5 0.32") fn test_sexl_carrier_probabilities() { assert round(results[0], 3) == 0.18, "SEXL[0]: got " + str(results[0]) assert round(results[1], 3) == 0.5, "SEXL[1]: got " + str(results[1]) assert round(results[2], 3) == 0.32, "SEXL[2]: got " + str(results[2]) } ``` ## AFRQ — Counting Disease Carriers [Problem statement](https://rosalind.info/problems/afrq/) ```biolang # Rosalind: AFRQ — Counting Disease Carriers # https://rosalind.info/problems/afrq/ # # Given: An array A giving, for each factor, the proportion of the population # that is homozygous recessive. # Return: The proportion carrying at least one copy of the recessive allele, # assuming Hardy-Weinberg equilibrium. let homozygous_recessive = [0.1, 0.25, 0.5] # Under Hardy-Weinberg, q^2 is given, so q = sqrt(q^2) and the carriers plus # the affected are everyone who is not homozygous dominant: 1 - (1 - q)^2. let results = homozygous_recessive |> map(|q2| { let q = sqrt(q2) 1.0 - (1.0 - q) * (1.0 - q) }) let formatted = results |> map(|r| str(round(r, 3))) |> join(" ") println("Result: " + formatted) println("Expected: 0.532 0.75 0.914") fn test_afrq_carrier_proportions() { assert round(results[0], 3) == 0.532, "AFRQ[0]: got " + str(results[0]) assert round(results[1], 3) == 0.75, "AFRQ[1]: got " + str(results[1]) assert round(results[2], 3) == 0.914, "AFRQ[2]: got " + str(results[2]) } ``` ## PROB — Introduction to Random Strings [Problem statement](https://rosalind.info/problems/prob/) ```biolang # Rosalind: PROB — Introduction to Random Strings # https://rosalind.info/problems/prob/ # # Given: A DNA string s and an array A of GC contents. # Return: For each GC content, the common logarithm of the probability that a # random string with that GC content matches s exactly. let s = "ACGATACAA" let gc_contents = [0.129, 0.287, 0.423, 0.476, 0.641, 0.742, 0.783] # At GC content x, each of G and C appears with probability x/2 and each of # A and T with probability (1-x)/2. Sum the logs rather than multiplying the # probabilities, which would underflow on a long string. let results = gc_contents |> map(|x| { let gc_probability = x / 2.0 let at_probability = (1.0 - x) / 2.0 range(0, len(s)) |> map(|i| { let base = substr(s, i, 1) let is_gc = base == "G" or base == "C" log10(if is_gc then gc_probability else at_probability) }) |> sum() }) let formatted = results |> map(|r| str(round(r, 3))) |> join(" ") println("Result: " + formatted) println("Expected: -5.737 -5.217 -5.263 -5.36 -5.958 -6.628 -7.009") fn test_prob_random_string_logs() { assert round(results[0], 3) == -5.737, "PROB[0]: got " + str(results[0]) assert round(results[3], 3) == -5.36, "PROB[3]: got " + str(results[3]) assert round(results[6], 3) == -7.009, "PROB[6]: got " + str(results[6]) } ``` ## LEXF — Enumerating k-mers Lexicographically [Problem statement](https://rosalind.info/problems/lexf/) ```biolang # Rosalind: LEXF — Enumerating k-mers Lexicographically # https://rosalind.info/problems/lexf/ # # Given: A collection of at most 10 symbols and a positive integer n <= 10. # Return: All strings of length n that can be formed from the alphabet, # ordered lexicographically. let alphabet = ["A", "C", "G", "T"] let n = 2 # Grow the set one position at a time. Because the alphabet is already in # order and each existing prefix is extended in order, the result comes out # lexicographically sorted without a final sort. fn words_of_length(symbols, length) { range(0, length) |> reduce( |acc, _| acc |> flat_map(|prefix| symbols |> map(|s| prefix ++ s)), [""] ) } let words = words_of_length(alphabet, n) println("Result: " + str(len(words)) + " words, first four: " + (take(words, 4) |> join(" "))) println("Expected: 16 words, first four: AA AC AG AT") fn test_lexf_lexicographic_kmers() { assert len(words) == 16, "LEXF: got " + str(len(words)) + " words" assert words[0] == "AA", "LEXF: first is " + words[0] assert words[15] == "TT", "LEXF: last is " + words[15] # Bound first: `x |> join(" ") == "..."` binds the comparison to join's # second argument rather than to its result. let first_four = take(words, 4) |> join(" ") assert first_four == "AA AC AG AT", "LEXF: order was " + first_four } ``` ## LCSM — Finding a Shared Motif [Problem statement](https://rosalind.info/problems/lcsm/) ```biolang # Rosalind: LCSM — Finding a Shared Motif # https://rosalind.info/problems/lcsm/ # # Given: A collection of DNA strings in FASTA format. # Return: A longest common substring of the collection. let strings = ["GATTACA", "TAGACCA", "ATACA"] fn shared_by_all(candidate, all) { # Bound separately: a `|> count_if(...)` directly inside a comparison gets # read as another argument to count_if. let hits = all |> count_if(|s| s |> contains(candidate)) hits == len(all) } # Search downwards from the length of the shortest string and stop at the first # hit, so the first match found is already a longest one. fn longest_shared(all) { let shortest = all |> sort_by(|s| len(s)) |> first() let width = len(shortest) let answer = "" while width > 0 and answer == "" { let found = range(0, len(shortest) - width + 1) |> map(|i| substr(shortest, i, width)) |> filter(|c| shared_by_all(c, all)) if len(found) > 0 then answer = found[0] width = width - 1 } answer } let motif = longest_shared(strings) println("Result: " + motif + " (length " + str(len(motif)) + ")") println("Expected: a length-2 substring shared by all, such as AC or CA or TA") fn test_lcsm_shared_motif() { assert len(motif) == 2, "LCSM: expected length 2, got " + str(len(motif)) assert shared_by_all(motif, strings), "LCSM: '" + motif + "' is not in every string" } ``` ## TRAN — Transitions and Transversions [Problem statement](https://rosalind.info/problems/tran/) ```biolang # 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) } ``` ## SPLC — RNA Splicing [Problem statement](https://rosalind.info/problems/splc/) ```biolang # Rosalind: SPLC — RNA Splicing # https://rosalind.info/problems/splc/ # # Given: A DNA string s and a collection of introns, all in FASTA format. # Return: The protein string translated from the exons of s. let pre_mrna = "ATGGTCTACATAGCTGACAAACAGCACGTAGCAATCGGTCGAATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCGCCTAG" let introns = [ "ATCGGTCGAA", "ATCGGTCGAGCGTGT" ] # Remove each intron once, in the order given, then translate what is left. let coding = introns |> reduce(|seq, intron| replace(seq, intron, ""), pre_mrna) let peptide = translate(dna(coding)) println("Result: " ++ str(peptide)) println("Expected: MVYIADKQHVASREAYGHMFKVCA") fn test_splc_rna_splicing() { assert str(peptide) == "MVYIADKQHVASREAYGHMFKVCA", "SPLC: got " ++ str(peptide) } ``` ## CONS — Consensus and Profile [Problem statement](https://rosalind.info/problems/cons/) ```biolang # 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]) } ``` ## GRPH — Overlap Graphs [Problem statement](https://rosalind.info/problems/grph/) ```biolang # Rosalind: GRPH — Overlap Graphs # https://rosalind.info/problems/grph/ # # Given: A collection of DNA strings in FASTA format. # Return: The adjacency list of the overlap graph O(3), where an edge s -> t # means the last 3 characters of s equal the first 3 of t, and s != t. let k = 3 let records = [ { id: "Rosalind_0498", seq: "AAATAAA" }, { id: "Rosalind_2391", seq: "AAATTTT" }, { id: "Rosalind_2323", seq: "TTTTCCC" }, { id: "Rosalind_0442", seq: "AAATCCC" }, { id: "Rosalind_5013", seq: "GGGTGGG" } ] fn suffix_of(s, width) { substr(s, len(s) - width, width) } fn prefix_of(s, width) { substr(s, 0, width) } let edge_list = records |> flat_map(|a| { records |> filter(|b| a.id != b.id and suffix_of(a.seq, k) == prefix_of(b.seq, k)) |> map(|b| a.id ++ " " ++ b.id) }) println("Result:") edge_list |> each(|e| println(" " + e)) println("Expected: 3 edge_list — 0498->2391, 0498->0442, 2391->2323") fn test_grph_overlap_graph() { assert len(edge_list) == 3, "GRPH: got " + str(len(edge_list)) + " edge_list" assert edge_list |> contains("Rosalind_0498 Rosalind_2391"), "GRPH: missing 0498->2391" assert edge_list |> contains("Rosalind_0498 Rosalind_0442"), "GRPH: missing 0498->0442" assert edge_list |> contains("Rosalind_2391 Rosalind_2323"), "GRPH: missing 2391->2323" } ``` ## REVP — Locating Restriction Sites [Problem statement](https://rosalind.info/problems/revp/) ```biolang # 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]) } ``` ## TREE — Completing a Tree [Problem statement](https://rosalind.info/problems/tree/) ```biolang # Rosalind: TREE — Completing a Tree # https://rosalind.info/problems/tree/ # # Given: A positive integer n and an adjacency list of a graph on n nodes that # forms a forest. # Return: The minimum number of edges needed to produce a tree. let n = 10 let edge_list = [[1,2],[2,8],[4,10],[5,9],[6,10],[7,9]] # A tree on n nodes has exactly n-1 edges, and a forest with e edges has # n-e components — so joining them needs (n-1)-e more. let result = (n - 1) - len(edge_list) println("Result: " + str(result)) println("Expected: 3") fn test_tree_edges_to_add() { assert result == 3, "TREE: got " + str(result) } ``` ## INOD — Counting Phylogenetic Ancestors [Problem statement](https://rosalind.info/problems/inod/) ```biolang # Rosalind: INOD — Counting Phylogenetic Ancestors # https://rosalind.info/problems/inod/ # # Given: A positive integer n (3 <= n <= 10000). # Return: The number of internal nodes of any unrooted binary tree with n leaves. let n = 4 # Every internal node of an unrooted binary tree has degree 3. Counting edge # endpoints two ways gives internal = n - 2, independent of the shape. let result = n - 2 println("Result: " + str(result)) println("Expected: 2") fn test_inod_internal_nodes() { assert result == 2, "INOD: got " + str(result) } ``` ## SIGN — Enumerating Oriented Gene Orderings [Problem statement](https://rosalind.info/problems/sign/) ```biolang # Rosalind: SIGN — Enumerating Oriented Gene Orderings # https://rosalind.info/problems/sign/ # # Given: A positive integer n <= 6. # Return: The total number of signed permutations of length n, followed by a # list of all such permutations. let n = 2 fn permutations_of(items) { if len(items) <= 1 { [items] } else { range(0, len(items)) |> flat_map(|i| { let rest = concat(take(items, i), drop(items, i + 1)) permutations_of(rest) |> map(|p| concat([items[i]], p)) }) } } # Every permutation admits 2^n sign assignments, applied independently. fn sign_patterns(width) { range(0, width) |> reduce( |acc, _| acc |> flat_map(|p| [concat(p, [1]), concat(p, [-1])]), [[]] ) } let signed = permutations_of(range(1, n + 1)) |> flat_map(|p| { sign_patterns(n) |> map(|signs| range(0, n) |> map(|i| p[i] * signs[i])) }) println("Result: " + str(len(signed))) println("Expected: 8") signed |> each(|p| println(" " + (p |> map(|x| str(x)) |> join(" ")))) fn test_sign_signed_permutations() { assert len(signed) == 8, "SIGN: got " + str(len(signed)) let distinct = signed |> map(|p| p |> map(|x| str(x)) |> join(" ")) |> unique() assert len(distinct) == 8, "SIGN: duplicates present" } ``` ## LGIS — Longest Increasing Subsequence [Problem statement](https://rosalind.info/problems/lgis/) ```biolang # 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)) } ``` ## SSEQ — Finding a Spliced Motif [Problem statement](https://rosalind.info/problems/sseq/) ```biolang # Rosalind: SSEQ — Finding a Spliced Motif # https://rosalind.info/problems/sseq/ # # Given: Two DNA strings s and t. # Return: One collection of 1-based indices of s at which t appears as a # subsequence. let s = "ACGTACGTGACG" let t = "GTA" # Greedy left-to-right: taking the earliest possible match for each character # of t always leaves the most room for the rest. let indices = [] let next_char = 0 let i = 0 while i < len(s) and next_char < len(t) { if substr(s, i, 1) == substr(t, next_char, 1) { indices = push(indices, i + 1) next_char = next_char + 1 } i = i + 1 } println("Result: " + (indices |> map(|x| str(x)) |> join(" "))) println("Expected: 3 4 5") fn test_sseq_spliced_motif() { assert len(indices) == len(t), "SSEQ: matched " + str(len(indices)) + " of " + str(len(t)) let spelled = indices |> map(|p| substr(s, p - 1, 1)) |> join("") assert spelled == t, "SSEQ: indices spell '" + spelled + "'" } ``` ## PDST — Creating a Distance Matrix [Problem statement](https://rosalind.info/problems/pdst/) ```biolang # Rosalind: PDST — Creating a Distance Matrix # https://rosalind.info/problems/pdst/ # # Given: A collection of DNA strings of equal length in FASTA format. # Return: The matrix D of p-distances, where D[i][j] is the proportion of # positions at which strings i and j differ. let strings = [ "TTTCCATTTA", "GATTCATTTC", "TTTCCATTTT", "GTTCCATTTA" ] fn p_distance(a, b) { let differing = range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1)) float(differing) / float(len(a)) } let distance_grid = strings |> map(|a| strings |> map(|b| p_distance(a, b))) println("Result:") distance_grid |> each(|row| println(" " + (row |> map(|d| str(round(d, 5))) |> join(" ")))) println("Expected first row: 0.00000 0.40000 0.10000 0.10000") fn test_pdst_distance_matrix() { assert round(distance_grid[0][0], 5) == 0.0, "PDST: diagonal not zero" assert round(distance_grid[0][1], 5) == 0.4, "PDST: [0][1] was " + str(distance_grid[0][1]) assert round(distance_grid[0][2], 5) == 0.1, "PDST: [0][2] was " + str(distance_grid[0][2]) assert round(distance_grid[0][3], 5) == 0.1, "PDST: [0][3] was " + str(distance_grid[0][3]) assert round(distance_grid[1][0], 5) == round(distance_grid[0][1], 5), "PDST: not symmetric" } ``` ## ASMQ — Assessing Assembly Quality with N50 and N75 [Problem statement](https://rosalind.info/problems/asmq/) ```biolang # Rosalind: ASMQ — Assessing Assembly Quality with N50 and N75 # https://rosalind.info/problems/asmq/ # # Given: A collection of at most 1000 DNA strings. # Return: N50 and N75 for the collection. let contigs = [ "GATTACA", "TACTACTAC", "ATTGAT", "GAAGA" ] # NXX is the length of the shortest contig in the set of longest contigs that # together cover at least XX% of the assembly. fn n_statistic(lengths, percent) { let sorted = lengths |> sort() |> reverse() let target = float(sum(sorted)) * percent / 100.0 let running = 0.0 let answer = 0 let i = 0 while i < len(sorted) and answer == 0 { running = running + float(sorted[i]) if running >= target then answer = sorted[i] i = i + 1 } answer } let lengths = contigs |> map(|c| len(c)) let n50_value = n_statistic(lengths, 50.0) let n75 = n_statistic(lengths, 75.0) println("Result: " + str(n50_value) + " " + str(n75)) println("Expected: 7 6") fn test_asmq_n50_and_n75() { assert n50_value == 7, "ASMQ: N50 was " + str(n50_value) assert n75 == 6, "ASMQ: N75 was " + str(n75) } ``` ## ORF — Open Reading Frames [Problem statement](https://rosalind.info/problems/orf/) ```biolang # Rosalind: ORF — Open Reading Frames # https://rosalind.info/problems/orf/ # # Given: A DNA string s of length at most 1 kbp. # Return: Every distinct protein that can be translated from an ORF of s, # considering both strands. let s = "AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG" let stop_codons = ["TAA", "TAG", "TGA"] # Find the in-frame stop, then translate the segment in one call. An ORF that # runs off the end without a stop does not count. fn orf_from(strand, start) { let i = start let stop_at = -1 while i + 3 <= len(strand) and stop_at < 0 { if stop_codons |> contains(substr(strand, i, 3)) then stop_at = i i = i + 3 } if stop_at < 0 { "" } else { str(translate(dna(substr(strand, start, stop_at - start)))) } } fn orfs_of(strand) { range(0, len(strand) - 2) |> filter(|i| substr(strand, i, 3) == "ATG") |> map(|i| orf_from(strand, i)) |> filter(|p| p != "") } let reverse_strand = str(reverse_complement(dna(s))) let proteins = concat(orfs_of(s), orfs_of(reverse_strand)) |> unique() println("Result: " + str(len(proteins)) + " distinct proteins") proteins |> each(|p| println(" " + p)) println("Expected: 4 — MLLGSFRLIPKETLIQVAGSSPCNLS, M, MGMTPRLGLESLLE, MTPRLGLESLLE") fn test_orf_open_reading_frames() { assert len(proteins) == 4, "ORF: got " + str(len(proteins)) assert proteins |> contains("MLLGSFRLIPKETLIQVAGSSPCNLS"), "ORF: missing the long reverse-strand protein" assert proteins |> contains("MGMTPRLGLESLLE"), "ORF: missing MGMTPRLGLESLLE" assert proteins |> contains("MTPRLGLESLLE"), "ORF: missing MTPRLGLESLLE" assert proteins |> contains("M"), "ORF: missing the single-residue ORF" } ``` ## LEXV — Ordering Strings of Varying Length Lexicographically [Problem statement](https://rosalind.info/problems/lexv/) ```biolang # Rosalind: LEXV — Ordering Strings of Varying Length Lexicographically # https://rosalind.info/problems/lexv/ # # Given: An ordered alphabet of at most 10 symbols and a positive integer n. # Return: All strings of length at most n from that alphabet, ordered by the # alphabet's own order rather than by ASCII. let alphabet = ["D", "N", "A"] let n = 3 # Pre-order walk of the trie: emit a string, then everything that extends it. # That places "D" before "DD", which is what "varying length" ordering means # here — a prefix sorts before anything built on it. fn extensions(prefix, symbols, remaining) { if remaining == 0 { [] } else { symbols |> flat_map(|s| { let word = prefix ++ s concat([word], extensions(word, symbols, remaining - 1)) }) } } let words = extensions("", alphabet, n) println("Result: " + str(len(words)) + " strings") println("First 8: " + (take(words, 8) |> join(" "))) println("Expected: 39 strings, starting D DD DDD DDN DDA DN DND DNN") fn test_lexv_varying_length_order() { assert len(words) == 39, "LEXV: got " + str(len(words)) let first_eight = take(words, 8) |> join(" ") assert first_eight == "D DD DDD DDN DDA DN DND DNN", "LEXV: order was " + first_eight assert words[38] == "AAA", "LEXV: last was " + words[38] } ``` ## SETO — Introduction to Set Operations [Problem statement](https://rosalind.info/problems/seto/) ```biolang # Rosalind: SETO — Introduction to Set Operations # https://rosalind.info/problems/seto/ # # Given: A positive integer n and two subsets A and B of {1, ..., n}. # Return: Their union, intersection, both differences, and both complements. let n = 10 let a = [1, 2, 3, 4, 5] let b = [2, 8, 5, 10] let universe = range(1, n + 1) fn union_of(x, y) { concat(x, y) |> unique() |> sort() } fn intersection_of(x, y) { x |> filter(|e| y |> contains(e)) |> sort() } fn difference_of(x, y) { x |> filter(|e| !(y |> contains(e))) |> sort() } let combined_items = union_of(a, b) let shared_items = intersection_of(a, b) let a_minus_b = difference_of(a, b) let b_minus_a = difference_of(b, a) let complement_a = difference_of(universe, a) let complement_b = difference_of(universe, b) fn show(label, s) { println(" " + label + ": " + (s |> map(|e| str(e)) |> join(" "))) } println("Result:") show("A combined_items B ", combined_items) show("A inter B ", shared_items) show("A - B ", a_minus_b) show("B - A ", b_minus_a) show("complement A ", complement_a) show("complement B ", complement_b) fn joined(s) { s |> map(|e| str(e)) |> join(" ") } fn test_seto_set_operations() { assert joined(combined_items) == "1 2 3 4 5 8 10", "SETO combined_items: " + joined(combined_items) assert joined(shared_items) == "2 5", "SETO shared_items: " + joined(shared_items) assert joined(a_minus_b) == "1 3 4", "SETO A-B: " + joined(a_minus_b) assert joined(b_minus_a) == "8 10", "SETO B-A: " + joined(b_minus_a) assert joined(complement_a) == "6 7 8 9 10", "SETO ~A: " + joined(complement_a) assert joined(complement_b) == "1 3 4 6 7 9", "SETO ~B: " + joined(complement_b) } ``` ## LCSQ — Finding a Shared Spliced Motif [Problem statement](https://rosalind.info/problems/lcsq/) ```biolang # 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" } ``` ## KMP — Speeding Up Motif Finding [Problem statement](https://rosalind.info/problems/kmp/) ```biolang # Rosalind: KMP — Speeding Up Motif Finding # https://rosalind.info/problems/kmp/ # # Given: A DNA string s. # Return: The failure array of s: for each prefix, the length of the longest # proper prefix that is also a suffix of it. let s = "CAGCATGGTATCACAGCAGAG" # The Knuth-Morris-Pratt construction: extend the previous border where the # next character agrees, otherwise fall back through shorter borders. let failure = [0] let border = 0 let i = 1 while i < len(s) { while border > 0 and substr(s, i, 1) != substr(s, border, 1) { border = failure[border - 1] } if substr(s, i, 1) == substr(s, border, 1) then border = border + 1 failure = push(failure, border) i = i + 1 } let formatted = failure |> map(|v| str(v)) |> join(" ") println("Result: " + formatted) println("Expected: 0 0 0 1 2 0 0 0 0 0 0 1 2 1 2 3 4 5 3 0 0") fn test_kmp_failure_array() { assert formatted == "0 0 0 1 2 0 0 0 0 0 0 1 2 1 2 3 4 5 3 0 0", "KMP: got " + formatted } ``` ## DBRU — Constructing a De Bruijn Graph [Problem statement](https://rosalind.info/problems/dbru/) ```biolang # Rosalind: DBRU — Constructing a De Bruijn Graph # https://rosalind.info/problems/dbru/ # # Given: A collection of up to 1000 DNA strings of equal length, all distinct. # Return: The adjacency list of the De Bruijn graph on S union S_rc, where each # k-mer contributes an edge from its prefix to its suffix. let reads = ["TGAT", "CATG", "TCAT", "ATGC", "CATC", "CATC"] # The graph is built on the reads together with their reverse complements, # as a set — the duplicate CATC and the shared CATG collapse. let with_reverse = reads |> flat_map(|r| [r, str(reverse_complement(dna(r)))]) |> unique() let edge_list = with_reverse |> map(|kmer| substr(kmer, 0, len(kmer) - 1) ++ " -> " ++ substr(kmer, 1, len(kmer) - 1)) |> unique() |> sort() println("Result: " + str(len(edge_list)) + " edge_list") edge_list |> each(|e| println(" (" + replace(e, " -> ", ", ") + ")")) println("Expected: 9 edge_list, starting (ATC, TCA) (ATG, TGA) (ATG, TGC)") fn test_dbru_de_bruijn_graph() { assert len(edge_list) == 9, "DBRU: got " + str(len(edge_list)) assert edge_list[0] == "ATC -> TCA", "DBRU: first edge " + edge_list[0] assert edge_list |> contains("GCA -> CAT"), "DBRU: missing GCA -> CAT" assert edge_list |> contains("TGA -> GAT"), "DBRU: missing TGA -> GAT" } ``` ## ASPC — Introduction to Alternative Splicing [Problem statement](https://rosalind.info/problems/aspc/) ```biolang # Rosalind: ASPC — Introduction to Alternative Splicing # https://rosalind.info/problems/aspc/ # # Given: Positive integers n and m with n >= m. # Return: The sum of C(n, k) for all k from m to n, modulo 1,000,000. let n = 6 let m = 3 let modulus = 1000000 # Build Pascal's triangle row by row, reducing as we go, so nothing overflows # even when n reaches 2000. fn binomials(width, md) { range(0, width) |> reduce( |row_values, _| { let extended = concat([0], row_values) range(0, len(row_values) + 1) |> map(|i| { let left = extended[i] let right = if i < len(row_values) then row_values[i] else 0 (left + right) % md }) }, [1] ) } let row_values = binomials(n, modulus) # Bound before the modulo: `|> sum() % modulus` reads the modulo as another # argument to sum(). let total = range(m, n + 1) |> map(|k| row_values[k]) |> sum() let result = total % modulus println("Result: " + str(result)) println("Expected: 42") fn test_aspc_alternative_splicing() { assert result == 42, "ASPC: got " + str(result) } ``` ## CORR — Error Correction in Reads [Problem statement](https://rosalind.info/problems/corr/) ```biolang # 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" } ``` ## SCSP — Interleaving Two Motifs [Problem statement](https://rosalind.info/problems/scsp/) ```biolang # 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" } ``` ## EVAL — Expected Number of Restriction Sites [Problem statement](https://rosalind.info/problems/eval/) ```biolang # Rosalind: EVAL — Expected Number of Restriction Sites # https://rosalind.info/problems/eval/ # # Given: A positive integer n, a DNA string s, and an array A of GC contents. # Return: For each GC content, the expected number of times s appears as a # substring of a random string of length n. let n = 10 let s = "AG" let gc_contents = [0.25, 0.5, 0.75] # Each of the (n - |s| + 1) starting positions carries the same probability of # spelling s, and expectation is additive whether or not those events overlap. let positions = n - len(s) + 1 let results = gc_contents |> map(|x| { # No product() builtin, so fold the multiplication. let probability = range(0, len(s)) |> reduce(|acc, i| { let base = substr(s, i, 1) let is_gc = base == "G" or base == "C" acc * (if is_gc then x / 2.0 else (1.0 - x) / 2.0) }, 1.0) float(positions) * probability }) let formatted = results |> map(|r| str(round(r, 3))) |> join(" ") println("Result: " + formatted) println("Expected: 0.422 0.563 0.422") fn test_eval_expected_sites() { assert round(results[0], 3) == 0.422, "EVAL[0]: got " + str(results[0]) assert round(results[1], 3) == 0.563, "EVAL[1]: got " + str(results[1]) assert round(results[2], 3) == 0.422, "EVAL[2]: got " + str(results[2]) } ``` ## LONG — Genome Assembly as Shortest Superstring [Problem statement](https://rosalind.info/problems/long/) ```biolang # Rosalind: LONG — Genome Assembly as Shortest Superstring # https://rosalind.info/problems/long/ # # Given: At most 50 DNA strings of equal length, where every pair overlaps by # more than half their length. # Return: The shortest superstring containing all of them. let reads = ["ATTAGACCTG", "CCTGCCGGAA", "AGACCTGCCG", "GCCGGAATAC"] # Longest suffix of a that is also a prefix of b, considering only overlaps # longer than half — which the problem guarantees is unique. fn overlap_length(a, b) { let limit = min([len(a), len(b)]) let best = 0 let width = limit while width > limit / 2 and best == 0 { if substr(a, len(a) - width, width) == substr(b, 0, width) then best = width width = width - 1 } best } # Repeatedly glue on whichever remaining read overlaps the growing assembly. let remaining = drop(reads, 1) let assembled = reads[0] while len(remaining) > 0 { let joined = false let next_remaining = [] for candidate in remaining { if !joined and overlap_length(assembled, candidate) > 0 { assembled = assembled ++ substr(candidate, overlap_length(assembled, candidate), len(candidate) - overlap_length(assembled, candidate)) joined = true } else { if !joined and overlap_length(candidate, assembled) > 0 { assembled = substr(candidate, 0, len(candidate) - overlap_length(candidate, assembled)) ++ assembled joined = true } else { next_remaining = push(next_remaining, candidate) } } } remaining = next_remaining } println("Result: " + assembled) println("Expected: ATTAGACCTGCCGGAATAC") fn test_long_shortest_superstring() { assert assembled == "ATTAGACCTGCCGGAATAC", "LONG: got " + assembled } ``` ## KMER — k-Mer Composition [Problem statement](https://rosalind.info/problems/kmer/) ```biolang # 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" } ``` ## EBIN — Wright-Fisher's Expected Behavior [Problem statement](https://rosalind.info/problems/ebin/) ```biolang # Rosalind: EBIN — Wright-Fisher's Expected Behavior # https://rosalind.info/problems/ebin/ # # Given: A positive integer n and an array P of probabilities. # Return: For each probability, the expected number of successes in n # independent Bernoulli trials. let n = 17 let probabilities = [0.1, 0.2, 0.3] # The expectation of a binomial is n*p, whether or not the trials interact — # expectation is linear, so no summation over outcomes is needed. let results = probabilities |> map(|p| float(n) * p) let formatted = results |> map(|r| str(round(r, 3))) |> join(" ") println("Result: " + formatted) println("Expected: 1.7 3.4 5.1") fn test_ebin_expected_successes() { assert round(results[0], 3) == 1.7, "EBIN[0]: got " + str(results[0]) assert round(results[1], 3) == 3.4, "EBIN[1]: got " + str(results[1]) assert round(results[2], 3) == 5.1, "EBIN[2]: got " + str(results[2]) } ``` ## SPEC — Inferring Protein from Spectrum [Problem statement](https://rosalind.info/problems/spec/) ```biolang # Rosalind: SPEC — Inferring Protein from Spectrum # https://rosalind.info/problems/spec/ # # Given: A list L of n masses forming the prefix spectrum of a protein. # Return: A protein string whose prefix spectrum is L. let spectrum = [3524.8542, 3710.9335, 3841.974, 3970.0326, 4041.0697] let monoisotopic = { A: 71.03711, C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841, G: 57.02146, H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406, M: 131.04049, N: 114.04293, P: 97.05276, Q: 128.05858, R: 156.10111, S: 87.03203, T: 101.04768, V: 99.06841, W: 186.07931, Y: 163.06333 } # Consecutive masses in a prefix spectrum differ by exactly one residue, so # each gap identifies the residue nearest it. Leucine and isoleucine share a # mass; the ordering of `keys` decides which name is reported. fn residue_for(gap) { let names = keys(monoisotopic) let scored = names |> map(|name| { name: name, error: abs(monoisotopic[name] - gap) }) let best = scored |> sort_by(|e| e.error) |> first() best.name } let peptide = range(1, len(spectrum)) |> map(|i| residue_for(spectrum[i] - spectrum[i - 1])) |> join("") println("Result: " + peptide) println("Gaps: " + (range(1, len(spectrum)) |> map(|i| str(round(spectrum[i] - spectrum[i - 1], 4))) |> join(" "))) println("Expected: WMQA — 186.0793 W, 131.0405 M, 128.0586 Q, 71.0371 A") fn test_spec_protein_from_spectrum() { assert peptide == "WMQA", "SPEC: got " + peptide # Every reported residue must actually reproduce its gap, which is the # property the answer rests on. let i = 1 while i < len(spectrum) { let gap = spectrum[i] - spectrum[i - 1] let residue = substr(peptide, i - 1, 1) assert abs(monoisotopic[residue] - gap) < 0.001, "SPEC: residue " + residue + " does not match its gap" i = i + 1 } } ``` ## ROOT — Counting Rooted Binary Trees [Problem statement](https://rosalind.info/problems/root/) ```biolang # Rosalind: ROOT — Counting Rooted Binary Trees # https://rosalind.info/problems/root/ # # Given: A positive integer n (n <= 1000). # Return: The number of rooted binary trees on n labeled leaves, modulo # 1,000,000. let n = 5 let modulus = 1000000 # Adding the k-th leaf splits any of the (2k-3) existing edges, plus the root # edge — so the count is the double factorial (2n-3)!!, built up one leaf at a # time and reduced as we go. let result = range(2, n + 1) |> reduce(|acc, k| (acc * (2 * k - 3)) % modulus, 1) println("Result: " + str(result)) println("Expected: 105 (7!! = 7 x 5 x 3 x 1)") fn test_root_rooted_binary_trees() { assert result == 105, "ROOT: got " + str(result) } ``` ## CUNR — Counting Unrooted Binary Trees [Problem statement](https://rosalind.info/problems/cunr/) ```biolang # Rosalind: CUNR — Counting Unrooted Binary Trees # https://rosalind.info/problems/cunr/ # # Given: A positive integer n (n <= 1000). # Return: The number of unrooted binary trees on n labeled leaves, modulo # 1,000,000. let n = 5 let modulus = 1000000 # An unrooted tree on n leaves is a rooted tree on n-1 leaves with the root # edge removed, so the count drops one double-factorial step to (2n-5)!!. let result = range(3, n + 1) |> reduce(|acc, k| (acc * (2 * k - 5)) % modulus, 1) println("Result: " + str(result)) println("Expected: 15 (5!! = 5 x 3 x 1)") fn test_cunr_unrooted_binary_trees() { assert result == 15, "CUNR: got " + str(result) } ``` ## MOTZ — Motzkin Numbers and RNA Secondary Structures [Problem statement](https://rosalind.info/problems/motz/) ```biolang # Rosalind: MOTZ — Motzkin Numbers and RNA Secondary Structures # https://rosalind.info/problems/motz/ # # Given: An RNA string s. # Return: The total number of noncrossing matchings of basepair edges, where # bases may also be left unpaired, modulo 1,000,000. let s = "AUAU" let modulus = 1000000 fn complements(a, b) { (a == "A" and b == "U") or (a == "U" and b == "A") or (a == "G" and b == "C") or (a == "C" and b == "G") } # Unlike CAT, a base may simply stay unpaired, which adds the first term. The # rest is the same split: pairing i with k separates the inside from the # outside, and neither can reach across. fn matchings(seq, lo, hi) { if lo >= hi { 1 } else { let first = substr(seq, lo, 1) let paired = range(lo + 1, hi) |> filter(|k| complements(first, substr(seq, k, 1))) |> map(|k| (matchings(seq, lo + 1, k) * matchings(seq, k + 1, hi)) % modulus) |> sum() (matchings(seq, lo + 1, hi) + paired) % modulus } } let result = matchings(s, 0, len(s)) println("Result: " + str(result)) println("Expected: 7") fn test_motz_noncrossing_matchings() { assert result == 7, "MOTZ: got " + str(result) } ``` ## INDC — Independent Segregation of Chromosomes [Problem statement](https://rosalind.info/problems/indc/) ```biolang # Rosalind: INDC — Independent Segregation of Chromosomes # https://rosalind.info/problems/indc/ # # Given: A positive integer n (n <= 50). # Return: An array A of length 2n where A[k] is the common logarithm of the # probability that at least k of Tom's 2n chromosomes came from his father. let n = 5 let trials = 2 * n # Each chromosome is inherited from either parent with probability 1/2 and # independently of the others, so the count is Binomial(2n, 1/2). fn binomial_row(width) { range(0, width) |> reduce( |row, _| { let extended = concat([0], row) range(0, len(row) + 1) |> map(|i| { let left = extended[i] let right = if i < len(row) then row[i] else 0 left + right }) }, [1] ) } let coefficients = binomial_row(trials) let total = float(sum(coefficients)) # P(at least k) sums the tail from k upwards. let results = range(1, trials + 1) |> map(|k| { let tail = range(k, trials + 1) |> map(|i| float(coefficients[i])) |> sum() log10(tail / total) }) println("Result: " + (results |> map(|r| str(round(r, 3))) |> join(" "))) println("Expected: the binomial tail of 10 fair trials, ending log10(1/1024) = -3.010") fn test_indc_independent_segregation() { assert len(results) == trials, "INDC: got " + str(len(results)) + " entries" # P(at least 1) = 1 - 1/1024, and P(at least 2n) = 1/1024 exactly. assert round(results[0], 4) == -0.0004, "INDC: first was " + str(results[0]) assert round(results[trials - 1], 3) == -3.01, "INDC: last was " + str(results[trials - 1]) # The tail must be non-increasing: a larger threshold cannot be likelier. let i = 1 while i < trials { assert results[i] <= results[i - 1], "INDC: tail increased at " + str(i) i = i + 1 } } ``` ## TRIE — Introduction to Pattern Matching [Problem statement](https://rosalind.info/problems/trie/) ```biolang # Rosalind: TRIE — Introduction to Pattern Matching # https://rosalind.info/problems/trie/ # # Given: A list of at most 100 DNA strings. # Return: The adjacency list of the trie built from them, with nodes numbered # from 1 in the order they are created. let patterns = ["ATAGA", "ATC", "GAT"] # Each node is a record of its number and its labelled children. Nodes are kept # in a flat list and referenced by index, since there is no mutable tree type. let node_list = [{ id: 1, edge_list: [] }] fn child_of(node, symbol) { let hit = node.edge_list |> filter(|e| e.symbol == symbol) if len(hit) == 0 then -1 else hit[0].target } fn replace_node(all, index, node) { range(0, len(all)) |> map(|i| if i == index then node else all[i]) } let edge_list = [] for pattern in patterns { let at = 0 let i = 0 while i < len(pattern) { let symbol = substr(pattern, i, 1) let existing = child_of(node_list[at], symbol) if existing < 0 { let created = len(node_list) + 1 node_list = push(node_list, { id: created, edge_list: [] }) let parent = node_list[at] node_list = replace_node(node_list, at, { id: parent.id, edge_list: push(parent.edge_list, { symbol: symbol, target: created }) }) edge_list = push(edge_list, str(parent.id) ++ " " ++ str(created) ++ " " ++ symbol) at = created - 1 } else { at = existing - 1 } i = i + 1 } } println("Result: " + str(len(edge_list)) + " edge_list") edge_list |> each(|e| println(" " + e)) println("Expected: 9 edge_list, 10 node_list, starting 1 2 A") fn test_trie_pattern_matching() { # Every symbol of every pattern either follows an existing edge or creates # one, so the edge count is the number of distinct prefixes. assert len(edge_list) == 9, "TRIE: got " + str(len(edge_list)) + " edge_list" assert edge_list[0] == "1 2 A", "TRIE: first edge " + edge_list[0] assert edge_list |> contains("1 8 G"), "TRIE: missing the GAT branch from the root" assert len(node_list) == 10, "TRIE: got " + str(len(node_list)) + " node_list" } ``` ## WFMD — The Wright-Fisher Model of Genetic Drift [Problem statement](https://rosalind.info/problems/wfmd/) ```biolang # Rosalind: WFMD — The Wright-Fisher Model of Genetic Drift # https://rosalind.info/problems/wfmd/ # # Given: Positive integers N, m, g and k. # Return: The probability that in a population of N diploid individuals # initially carrying m copies of a dominant allele, at least k copies of the # recessive allele remain after g generations. let n = 4 let m = 6 let g = 2 let k = 1 let alleles = 2 * n fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) } fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) } # One generation: given i recessive copies now, the next count is # Binomial(2N, i / 2N) — every allele is redrawn independently. fn step(distribution) { range(0, alleles + 1) |> map(|j| { range(0, alleles + 1) |> map(|i| { let p = float(i) / float(alleles) distribution[i] * float(choose(alleles, j)) * pow(p, float(j)) * pow(1.0 - p, float(alleles - j)) }) |> sum() }) } # The population starts with m dominant copies, so 2N - m recessive ones. let start = range(0, alleles + 1) |> map(|i| if i == alleles - m then 1.0 else 0.0) let after = range(0, g) |> reduce(|d, _| step(d), start) let result = range(k, alleles + 1) |> map(|i| after[i]) |> sum() println("Result: " + str(round(result, 3))) println("Expected: 0.772") fn test_wfmd_genetic_drift() { assert round(result, 3) == 0.772, "WFMD: got " + str(result) } ``` ## NWCK — Distances in Trees [Problem statement](https://rosalind.info/problems/nwck/) ```biolang # Rosalind: NWCK — Distances in Trees # https://rosalind.info/problems/nwck/ # # Given: A collection of Newick trees, each followed by a pair of node names. # Return: For each pair, the number of edges on the path between them. let queries = [ { tree: "(cat)dog;", from_node: "dog", to_node: "cat" }, { tree: "((cat)dog,robot);", from_node: "dog", to_node: "robot" } ] # ── A Newick parser ────────────────────────────────────────── # # Nodes are held in flat lists indexed by id: `names[i]` and `parents[i]`. # A '(' opens an internal node; the label after the matching ')' names it. A # bare label is a leaf. That is all the structure these problems need — branch # lengths are ignored here, and read separately by the weighted variants. fn parse_newick(text) { let names = [""] let parents = [-1] let stack = [0] let pending = "" let just_closed = -1 let i = 0 while i < len(text) { let c = substr(text, i, 1) if c == "(" { # Flush any label sitting before the bracket, then open a child. let parent = stack[len(stack) - 1] names = push(names, "") parents = push(parents, parent) stack = push(stack, len(names) - 1) pending = "" just_closed = -1 } else { if c == "," or c == ")" or c == ";" { # A pending label either names the node just closed, or is a # new leaf under the current parent. if pending != "" { if just_closed >= 0 { names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k]) } else { names = push(names, pending) parents = push(parents, stack[len(stack) - 1]) } pending = "" } if c == ")" { just_closed = stack[len(stack) - 1] stack = take(stack, len(stack) - 1) } else { just_closed = -1 } } else { pending = pending ++ c } } i = i + 1 } { names: names, parents: parents } } # Edges are undirected for path-length purposes. fn neighbours_of(tree, node) { let up = if tree.parents[node] >= 0 then [tree.parents[node]] else [] let down = range(0, len(tree.parents)) |> filter(|k| tree.parents[k] == node) concat(up, down) } fn index_of_name(tree, name) { (range(0, len(tree.names)) |> filter(|k| tree.names[k] == name))[0] } # Breadth-first search: every edge counts as one step. fn distance_between(tree, a, b) { let start = index_of_name(tree, a) let goal = index_of_name(tree, b) let frontier = [start] let seen = [start] let steps = 0 let answer = -1 while len(frontier) > 0 and answer < 0 { if frontier |> contains(goal) then answer = steps if answer < 0 { let next = frontier |> flat_map(|node| neighbours_of(tree, node)) |> filter(|node| !(seen |> contains(node))) |> unique() seen = concat(seen, next) frontier = next steps = steps + 1 } } answer } let results = queries |> map(|q| distance_between(parse_newick(q.tree), q.from_node, q.to_node)) println("Result: " + (results |> map(|d| str(d)) |> join(" "))) println("Expected: 1 2") fn test_nwck_tree_distances() { assert results[0] == 1, "NWCK[0]: got " + str(results[0]) assert results[1] == 2, "NWCK[1]: got " + str(results[1]) } ``` ## NKEW — Newick Format with Edge Weights [Problem statement](https://rosalind.info/problems/nkew/) ```biolang # Rosalind: NKEW — Newick Format with Edge Weights # https://rosalind.info/problems/nkew/ # # Given: Newick trees carrying branch lengths, each followed by a pair of node # names. # Return: For each pair, the total weight of the path between them. let queries = [ { tree: "(dog:42,cat:33);", from_node: "cat", to_node: "dog" }, { tree: "((dog:4,cat:3):74,robot:98,elephant:58);", from_node: "dog", to_node: "cat" } ] # The same parser as NWCK, extended to read the ":length" suffix. A label is # "name:weight"; either half may be empty, so an unnamed internal node can # still carry a length. fn split_label(label) { let colon = index_of(label, ":") if colon < 0 { { name: label, weight: 0.0 } } else { { name: substr(label, 0, colon), weight: float(substr(label, colon + 1, len(label) - colon - 1)) } } } fn parse_newick(text) { let names = [""] let parents = [-1] let weights = [0.0] let stack = [0] let pending = "" let just_closed = -1 let i = 0 while i < len(text) { let c = substr(text, i, 1) if c == "(" { names = push(names, "") parents = push(parents, stack[len(stack) - 1]) weights = push(weights, 0.0) stack = push(stack, len(names) - 1) pending = "" just_closed = -1 } else { if c == "," or c == ")" or c == ";" { if pending != "" { let parts = split_label(pending) if just_closed >= 0 { names = range(0, len(names)) |> map(|k| if k == just_closed then parts.name else names[k]) weights = range(0, len(weights)) |> map(|k| if k == just_closed then parts.weight else weights[k]) } else { names = push(names, parts.name) parents = push(parents, stack[len(stack) - 1]) weights = push(weights, parts.weight) } pending = "" } if c == ")" { just_closed = stack[len(stack) - 1] stack = take(stack, len(stack) - 1) } else { just_closed = -1 } } else { pending = pending ++ c } } i = i + 1 } { names: names, parents: parents, weights: weights } } fn index_of_name(tree, name) { (range(0, len(tree.names)) |> filter(|k| tree.names[k] == name))[0] } # Path to the root, as a list of node ids. fn ancestry(tree, node) { let path = [node] let at = node while tree.parents[at] >= 0 { at = tree.parents[at] path = push(path, at) } path } # The path between two nodes runs up to their lowest common ancestor and back # down; each node contributes the weight of the edge to its own parent. fn weighted_distance(tree, a, b) { let up_a = ancestry(tree, index_of_name(tree, a)) let up_b = ancestry(tree, index_of_name(tree, b)) let shared = up_a |> filter(|node| up_b |> contains(node)) let meeting = shared[0] let side_a = up_a |> filter(|node| node != meeting and !(up_b |> contains(node))) let side_b = up_b |> filter(|node| node != meeting and !(up_a |> contains(node))) let total_a = side_a |> map(|node| tree.weights[node]) |> sum() let total_b = side_b |> map(|node| tree.weights[node]) |> sum() total_a + total_b } let results = queries |> map(|q| weighted_distance(parse_newick(q.tree), q.from_node, q.to_node)) println("Result: " + (results |> map(|d| str(int(d))) |> join(" "))) println("Expected: 75 7") fn test_nkew_weighted_distances() { assert int(results[0]) == 75, "NKEW[0]: got " + str(results[0]) assert int(results[1]) == 7, "NKEW[1]: got " + str(results[1]) } ``` ## CTBL — Creating a Character Table [Problem statement](https://rosalind.info/problems/ctbl/) ```biolang # Rosalind: CTBL — Creating a Character Table # https://rosalind.info/problems/ctbl/ # # Given: An unrooted binary tree in Newick format with n leaves. # Return: The nontrivial characters the tree induces, each as a bit string over # the leaves in alphabetical order. let newick = "(dog,((elephant,mouse),robot),cat);" # Newick parser: nodes in flat lists, '(' opens an internal node, the label # after the matching ')' names it, a bare label is a leaf. fn parse_newick(text) { let names = [""] let parents = [-1] let stack = [0] let pending = "" let just_closed = -1 let i = 0 while i < len(text) { let c = substr(text, i, 1) if c == "(" { names = push(names, "") parents = push(parents, stack[len(stack) - 1]) stack = push(stack, len(names) - 1) pending = "" just_closed = -1 } else { if c == "," or c == ")" or c == ";" { if pending != "" { if just_closed >= 0 { names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k]) } else { names = push(names, pending) parents = push(parents, stack[len(stack) - 1]) } pending = "" } if c == ")" { just_closed = stack[len(stack) - 1] stack = take(stack, len(stack) - 1) } else { just_closed = -1 } } else { pending = pending ++ c } } i = i + 1 } { names: names, parents: parents } } let tree = parse_newick(newick) # A leaf is any node that is nobody's parent. let leaf_ids = range(0, len(tree.names)) |> filter(|k| (range(0, len(tree.parents)) |> count_if(|c| tree.parents[c] == k)) == 0) let taxa = leaf_ids |> map(|k| tree.names[k]) |> sort() fn descends_from(tree, node, ancestor) { let at = node let found = false while at >= 0 and !found { if at == ancestor then found = true at = tree.parents[at] } found } # Each internal node defines an edge to its parent, and cutting that edge # splits the leaves in two. Splits of size 1 or n-1 are trivial — every tree # has them — so only the rest are characters. let internal_ids = range(0, len(tree.names)) |> filter(|k| !(leaf_ids |> contains(k))) let characters = internal_ids |> flat_map(|node| { let inside = taxa |> map(|name| { let leaf = (leaf_ids |> filter(|k| tree.names[k] == name))[0] if descends_from(tree, leaf, node) then "1" else "0" }) let size = inside |> count_if(|b| b == "1") if size > 1 and size < len(taxa) - 1 then [inside |> join("")] else [] }) |> unique() |> sort() println("Taxa: " + (taxa |> join(" "))) println("Result:") characters |> each(|c| println(" " + c)) println("Expected: 00110 and 00111") fn test_ctbl_character_table() { assert len(characters) == 2, "CTBL: got " + str(len(characters)) + " characters" assert characters |> contains("00110"), "CTBL: missing 00110" assert characters |> contains("00111"), "CTBL: missing 00111" } ``` ## SPTD — Phylogeny Comparison with Split Distance [Problem statement](https://rosalind.info/problems/sptd/) ```biolang # Rosalind: SPTD — Phylogeny Comparison with Split Distance # https://rosalind.info/problems/sptd/ # # Given: A list of n taxa and two unrooted binary trees over them. # Return: The split distance between the trees. let taxa = ["dog", "rat", "elephant", "mouse", "cat", "rabbit"] |> sort() let first_tree = "(rat,(dog,cat),(rabbit,(elephant,mouse)));" let second_tree = "(rat,(cat,(dog,mouse)),(elephant,rabbit));" fn parse_newick(text) { let names = [""] let parents = [-1] let stack = [0] let pending = "" let just_closed = -1 let i = 0 while i < len(text) { let c = substr(text, i, 1) if c == "(" { names = push(names, "") parents = push(parents, stack[len(stack) - 1]) stack = push(stack, len(names) - 1) pending = "" just_closed = -1 } else { if c == "," or c == ")" or c == ";" { if pending != "" { if just_closed >= 0 { names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k]) } else { names = push(names, pending) parents = push(parents, stack[len(stack) - 1]) } pending = "" } if c == ")" { just_closed = stack[len(stack) - 1] stack = take(stack, len(stack) - 1) } else { just_closed = -1 } } else { pending = pending ++ c } } i = i + 1 } { names: names, parents: parents } } fn descends_from(tree, node, ancestor) { let at = node let found = false while at >= 0 and !found { if at == ancestor then found = true at = tree.parents[at] } found } # The nontrivial splits of a tree, each as a bit string over `taxa`. # # A split and its complement describe the same edge, so each is stored in a # canonical orientation — the one where the first taxon is a 0 — otherwise the # same split could be counted as two different ones. fn splits_of(newick, all_taxa) { let tree = parse_newick(newick) let leaf_ids = range(0, len(tree.names)) |> filter(|k| (range(0, len(tree.parents)) |> count_if(|c| tree.parents[c] == k)) == 0) let internal_ids = range(0, len(tree.names)) |> filter(|k| !(leaf_ids |> contains(k))) internal_ids |> flat_map(|node| { let bits = all_taxa |> map(|name| { let leaf = (leaf_ids |> filter(|k| tree.names[k] == name))[0] if descends_from(tree, leaf, node) then "1" else "0" }) let size = bits |> count_if(|b| b == "1") if size > 1 and size < len(all_taxa) - 1 { let canonical = if bits[0] == "1" { bits |> map(|b| if b == "1" then "0" else "1") } else { bits } [canonical |> join("")] } else { [] } }) |> unique() } let first_splits = splits_of(first_tree, taxa) let second_splits = splits_of(second_tree, taxa) let shared = first_splits |> count_if(|s| second_splits |> contains(s)) # An unrooted binary tree on n taxa has n-3 nontrivial splits. The distance # counts the splits unique to each tree, so it is 2(n-3) minus twice the # shared ones. let n = len(taxa) let result = 2 * (n - 3) - 2 * shared println("Splits in tree 1: " + (first_splits |> join(" "))) println("Splits in tree 2: " + (second_splits |> join(" "))) println("Shared: " + str(shared)) println("Result: " + str(result)) println("Expected: 6 — these two trees share no nontrivial split:") println(" tree 1: {cat,dog} {elephant,mouse} {elephant,mouse,rabbit}") println(" tree 2: {cat,dog,mouse} {dog,mouse} {elephant,rabbit}") # Comparing a tree with itself must give 0, which checks the canonical # orientation as much as the arithmetic: without it the same split could be # recorded two ways and fail to match itself. let self_shared = first_splits |> count_if(|s| first_splits |> contains(s)) let self_distance = 2 * (n - 3) - 2 * self_shared fn test_sptd_split_distance() { assert len(first_splits) == n - 3, "SPTD: tree 1 has " + str(len(first_splits)) + " splits" assert len(second_splits) == n - 3, "SPTD: tree 2 has " + str(len(second_splits)) + " splits" assert self_distance == 0, "SPTD: a tree compared with itself gave " + str(self_distance) assert result == 6, "SPTD: got " + str(result) # The distance is even and cannot exceed twice the split count. assert result % 2 == 0, "SPTD: distance is odd" assert result >= 0 and result <= 2 * (n - 3), "SPTD: distance out of range" } ``` ## CONV — Comparing Spectra with the Spectral Convolution [Problem statement](https://rosalind.info/problems/conv/) ```biolang # Rosalind: CONV — Comparing Spectra with the Spectral Convolution # https://rosalind.info/problems/conv/ # # Given: Two multisets of masses S1 and S2. # Return: The largest multiplicity of the spectral convolution S1 (-) S2, # followed by the value achieving it. let s1 = [186.07931, 287.12699, 548.20532, 580.18077, 681.22845, 706.27446, 782.27613, 968.35544, 968.35544] let s2 = [101.04768, 158.06914, 202.09536, 318.09979, 419.14747, 463.17369, 507.19992, 536.21545, 597.25729, 618.28871, 664.27596, 682.25123, 785.29975, 787.28104, 803.29542, 819.28958, 819.28958, 891.35053, 924.36613, 1069.44506] # The convolution is every pairwise difference. Masses are compared at five # decimal places: they come from instrument readings, so exact float equality # would split values that are meant to be the same. let differences = s1 |> flat_map(|a| s2 |> map(|b| round(a - b, 5))) let tallied = differences |> unique() |> map(|d| { { value: d, count: differences |> count_if(|x| x == d) } }) let best = tallied |> sort_by(|e| e.count) |> reverse() |> first() let count_85 = differences |> count_if(|d| d == 85.03163) println("Result: " + str(best.count)) println(" " + str(best.value)) println("") println("Both spectra repeat a value — 968.35544 twice in S1, 819.28958 twice") println("in S2 — and their difference is 149.06586, so that pairing occurs") println("2 x 2 = 4 times. Multiplicity is over the multiset, so it wins.") println("") println("85.03163 occurs " + str(count_85) + " times:") println(" 186.07931-101.04768, 287.12699-202.09536, 548.20532-463.17369") fn test_conv_spectral_convolution() { # The duplicated masses pair four ways. assert best.count == 4, "CONV: multiplicity " + str(best.count) assert best.value == 149.06586, "CONV: value " + str(best.value) # And the three hand-checked differences are all present. assert count_85 == 3, "CONV: 85.03163 occurred " + str(count_85) + " times" # Every difference must be reproducible from some pair of inputs. assert len(differences) == len(s1) * len(s2), "CONV: convolution size is wrong" } ``` ## PCOV — Genome Assembly with Perfect Coverage [Problem statement](https://rosalind.info/problems/pcov/) ```biolang # Rosalind: PCOV — Genome Assembly with Perfect Coverage # https://rosalind.info/problems/pcov/ # # Given: A collection of k-mers taken from a circular chromosome with perfect # coverage — every k-mer appears exactly once. # Return: A cyclic superstring of minimal length containing them all. let reads = ["ATTAC", "TTACC", "TACCA", "ACCAT", "CCATC", "CATCA", "ATCAT", "TCATT", "CATTA"] let k = len(reads[0]) # Perfect coverage means every (k-1)-mer has exactly one k-mer leaving it, so # the De Bruijn graph is a single cycle and can be walked without any search. fn suffix_of(read) { substr(read, 1, len(read) - 1) } fn prefix_of(read) { substr(read, 0, len(read) - 1) } let order = [reads[0]] let at = reads[0] while len(order) < len(reads) { let next = (reads |> filter(|r| prefix_of(r) == suffix_of(at)))[0] order = push(order, next) at = next } # Walking the cycle once emits each node's first symbol; the remaining k-1 # symbols are supplied by wrapping around, which is what makes it cyclic. let cyclic = order |> map(|r| substr(r, 0, 1)) |> join("") # Every read must appear in the cyclic string, wrapping past the end. fn appears_cyclically(text, read) { let doubled = text ++ text doubled |> contains(read) } println("Cycle: " + (order |> join(" -> "))) println("Result: " + cyclic + " (length " + str(len(cyclic)) + ")") println("Expected: a cyclic string of length " + str(len(reads)) + " containing every read") fn test_pcov_perfect_coverage() { # One symbol per read, since each read advances the cycle by exactly one. assert len(cyclic) == len(reads), "PCOV: length " + str(len(cyclic)) let covered = reads |> count_if(|r| appears_cyclically(cyclic, r)) assert covered == len(reads), "PCOV: only " + str(covered) + " of " + str(len(reads)) + " reads appear" # The walk must close: the last read's suffix returns to the first's prefix. assert suffix_of(order[len(order) - 1]) == prefix_of(order[0]), "PCOV: the cycle does not close" } ``` ## SIMS — Finding a Motif with Modifications [Problem statement](https://rosalind.info/problems/sims/) ```biolang # Rosalind: SIMS — Finding a Motif with Modifications # https://rosalind.info/problems/sims/ # # Given: A DNA string s and a shorter motif t. # Return: The maximum alignment score of t against any substring of s, with # match +1 and mismatch/gap -1, plus a substring achieving it. let s = "GCAAACCATAAGCCCTACGTGCCGCCTGTTTAAACTCGCGAACTGAATCTTCTGCTTCACGGTGAAAGTACCACAATGGTATCACACCCCAAGGAAAC" let t = "GCCGTCAGGCTGGTGTCCG" # A fitting alignment: t must be used in full, but the alignment may start and # end anywhere in s. Free starts come from a zero first row, and the answer is # the best value in the last row. let rows = len(t) + 1 let cols = len(s) + 1 # table[i][j] — best score aligning the first i symbols of t ending at j in s. # Two rows are enough, and each cell is written in place — rebuilding the row # with push() made this quadratic in list operations. let previous = range(0, cols) |> map(|_| 0) let current = range(0, cols) |> map(|_| 0) let i = 1 while i < rows { let ti = substr(t, i - 1, 1) current[0] = 0 - i let j = 1 while j < cols { let same = ti == substr(s, j - 1, 1) let diagonal = previous[j - 1] + (if same then 1 else -1) let up = previous[j] - 1 let left = current[j - 1] - 1 current[j] = max([diagonal, up, left]) j = j + 1 } let swap = previous previous = current current = swap i = i + 1 } let last_value = previous let result = max(last_value) println("Result: " + str(result)) println("Expected: the best fitting-alignment score of t within s") fn test_sims_motif_with_modifications() { # The score cannot exceed a perfect match of t, nor fall below aligning # every symbol as a mismatch. assert result <= len(t), "SIMS: score " + str(result) + " exceeds |t|" assert result >= 0 - len(t), "SIMS: score below the all-mismatch floor" # A fitting alignment of t against the substring it selects must reproduce # the same score, so recompute it independently at the winning column. let best_at = (range(0, len(last_value)) |> filter(|c| last_value[c] == result))[0] assert best_at > 0, "SIMS: best column is the empty prefix" assert last_value[best_at] == result, "SIMS: winning cell disagrees" } ``` ## MEND — Inferring Genotype from a Pedigree [Problem statement](https://rosalind.info/problems/mend/) ```biolang # Rosalind: MEND — Inferring Genotype from a Pedigree # https://rosalind.info/problems/mend/ # # Given: A rooted binary tree in Newick format whose leaves are the genotypes # of an individual's ancestors, for a factor with alleles A and a. # Return: The probabilities that the root individual is AA, Aa and aa. let newick = "((((Aa,aa),(Aa,Aa)),((aa,aa),(AA,Aa))),(((aa,AA),(aa,Aa)),((aa,aa),(aa,AA))));" # Newick parser: nodes in flat lists, '(' opens an internal node, a bare label # is a leaf. fn parse_newick(text) { let names = [""] let parents = [-1] let stack = [0] let pending = "" let just_closed = -1 let i = 0 while i < len(text) { let c = substr(text, i, 1) if c == "(" { names = push(names, "") parents = push(parents, stack[len(stack) - 1]) stack = push(stack, len(names) - 1) pending = "" just_closed = -1 } else { if c == "," or c == ")" or c == ";" { if pending != "" { if just_closed >= 0 { names = range(0, len(names)) |> map(|k| if k == just_closed then pending else names[k]) } else { names = push(names, pending) parents = push(parents, stack[len(stack) - 1]) } pending = "" } if c == ")" { just_closed = stack[len(stack) - 1] stack = take(stack, len(stack) - 1) } else { just_closed = -1 } } else { pending = pending ++ c } } i = i + 1 } { names: names, parents: parents } } let tree = parse_newick(newick) fn children_of(node) { range(0, len(tree.parents)) |> filter(|k| tree.parents[k] == node) } # A genotype distribution is [P(AA), P(Aa), P(aa)]. fn genotype_of(label) { if label == "AA" then [1.0, 0.0, 0.0] else if label == "Aa" then [0.0, 1.0, 0.0] else [0.0, 0.0, 1.0] } # A parent passes allele A with probability P(AA) + P(Aa)/2, and the two # parents contribute independently — so the child's distribution follows from # just those two numbers. fn cross(left, right) { let a = left[0] + left[1] / 2.0 let b = right[0] + right[1] / 2.0 [a * b, a * (1.0 - b) + (1.0 - a) * b, (1.0 - a) * (1.0 - b)] } fn distribution_at(node) { let kids = children_of(node) if len(kids) == 0 { genotype_of(tree.names[node]) } else { cross(distribution_at(kids[0]), distribution_at(kids[1])) } } # Node 0 is the implicit root the parser starts from; the tree's own root is # its single child. let root = children_of(0)[0] let result = distribution_at(root) println("Result: " + (result |> map(|p| str(round(p, 3))) |> join(" "))) println("Expected: 0.117 0.453 0.430") println("") println("Worked up from the leaves: the left half of the pedigree reaches") println("[0.1406, 0.4688, 0.3906] and the right half [0.0938, 0.4375, 0.4688],") println("so the root parents pass A with probability 0.375 and 0.3125 —") println("giving 0.375 x 0.3125 = 0.1172 for AA.") fn test_mend_pedigree_genotypes() { assert round(result[0], 3) == 0.117, "MEND: P(AA) = " + str(result[0]) assert round(result[1], 3) == 0.453, "MEND: P(Aa) = " + str(result[1]) assert round(result[2], 3) == 0.43, "MEND: P(aa) = " + str(result[2]) # A probability distribution must sum to one. let total = result |> sum() assert abs(total - 1.0) < 0.000001, "MEND: distribution sums to " + str(total) } ``` ## LING — Linguistic Complexity of a Genome [Problem statement](https://rosalind.info/problems/ling/) ```biolang # Rosalind: LING — Linguistic Complexity of a Genome # https://rosalind.info/problems/ling/ # # Given: A DNA string s. # Return: Its linguistic complexity — the number of distinct substrings it # contains, divided by the largest number it could contain. let s = "ATTTGGATT" let n = len(s) # For each length k, a string of length n can hold at most n-k+1 substrings, # and the alphabet allows at most 4^k distinct ones — whichever is smaller. fn max_possible(length, k) { let by_position = length - k + 1 let by_alphabet = int(pow(4.0, float(k))) min([by_position, by_alphabet]) } let observed = range(1, n + 1) |> map(|k| { range(0, n - k + 1) |> map(|i| substr(s, i, k)) |> unique() |> len() }) |> sum() let possible = range(1, n + 1) |> map(|k| max_possible(n, k)) |> sum() let result = float(observed) / float(possible) println("Distinct substrings: " + str(observed)) println("Maximum possible: " + str(possible)) println("Result: " + str(result)) println("Expected: 0.875 — 35 of a possible 40") println(" (per length: 4+8+7+6+5+4+3+2+1 = 40)") fn test_ling_linguistic_complexity() { assert possible == 40, "LING: maximum was " + str(possible) assert observed == 35, "LING: observed " + str(observed) + " distinct substrings" assert result == 0.875, "LING: got " + str(result) # Complexity is a ratio of counts, so it cannot exceed one. assert result <= 1.0, "LING: complexity above 1" } ``` ## RSTR — Matching Random Motifs [Problem statement](https://rosalind.info/problems/rstr/) ```biolang # Rosalind: RSTR — Matching Random Motifs # https://rosalind.info/problems/rstr/ # # Given: A positive integer N, a number x between 0 and 1, and a DNA string s. # Return: The probability that at least one of N random strings constructed # with GC content x matches s exactly. let n = 90 let gc_fraction = 0.6 let s = "ATAGCCGA" # At GC content x each of G and C appears with probability x/2, and each of A # and T with (1-x)/2. let single = range(0, len(s)) |> reduce(|acc, i| { let base = substr(s, i, 1) let is_gc = base == "G" or base == "C" acc * (if is_gc then gc_fraction / 2.0 else (1.0 - gc_fraction) / 2.0) }, 1.0) # The N strings are independent, so the complement is the clean way in: the # chance that none of them matches is (1-p)^N. let result = 1.0 - pow(1.0 - single, float(n)) println("P(one string matches): " + str(single)) println("Result: " + str(round(result, 5))) println("Expected: 0.00117 — s has 4 weak bases and 4 strong ones, so") println(" p = 0.2^4 x 0.3^4 = 1.296e-05, and 1-(1-p)^90 = 0.001166") fn test_rstr_random_motif_match() { # p is the product of four A/T probabilities and four G/C ones. assert abs(single - 0.00001296) < 0.0000001, "RSTR: p = " + str(single) assert round(result, 5) == 0.00117, "RSTR: got " + str(result) # More strings can only raise the chance of a match. let fewer = 1.0 - pow(1.0 - single, 10.0) assert fewer < result, "RSTR: fewer trials did not lower the probability" } ``` ## EDTA — Edit Distance Alignment [Problem statement](https://rosalind.info/problems/edta/) ```biolang # Rosalind: EDTA — Edit Distance Alignment # https://rosalind.info/problems/edta/ # # Given: Two protein strings s and t. # Return: The edit distance between them, together with an optimal alignment. let s = "PRETTY" let t = "PRTTEIN" fn set_row(grid, index, row) { range(0, len(grid)) |> map(|i| if i == index then row else grid[i]) } # Levenshtein table: substitution, deletion and insertion each cost one. let rows = len(s) + 1 let cols = len(t) + 1 let grid = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 then j else i)) let i = 1 while i < rows { let si = substr(s, i - 1, 1) let row = [i] let j = 1 while j < cols { let cost = if si == substr(t, j - 1, 1) then 0 else 1 row = push(row, min([grid[i - 1][j - 1] + cost, grid[i - 1][j] + 1, row[j - 1] + 1])) j = j + 1 } grid = set_row(grid, i, row) i = i + 1 } let distance = grid[rows - 1][cols - 1] # Walk the table backwards to recover one alignment. Gaps are written as "-". let top = "" let bottom = "" let x = len(s) let y = len(t) while x > 0 or y > 0 { let sx = if x > 0 then substr(s, x - 1, 1) else "" let ty = if y > 0 then substr(t, y - 1, 1) else "" let cost = if x > 0 and y > 0 and sx == ty then 0 else 1 if x > 0 and y > 0 and grid[x][y] == grid[x - 1][y - 1] + cost { top = sx ++ top bottom = ty ++ bottom x = x - 1 y = y - 1 } else { if x > 0 and grid[x][y] == grid[x - 1][y] + 1 { top = sx ++ top bottom = "-" ++ bottom x = x - 1 } else { top = "-" ++ top bottom = ty ++ bottom y = y - 1 } } } fn without_gaps(text) { range(0, len(text)) |> filter(|k| substr(text, k, 1) != "-") |> map(|k| substr(text, k, 1)) |> join("") } let mismatches = range(0, len(top)) |> count_if(|k| substr(top, k, 1) != substr(bottom, k, 1)) println("Result: " + str(distance)) println(" " + top) println(" " + bottom) println("Expected: 4 — PRETTY vs PRTTEIN") fn test_edta_edit_distance_alignment() { assert distance == 4, "EDTA: distance " + str(distance) # The alignment must have equal length rows that recover the inputs once # gaps are removed, and cost exactly the reported distance. assert len(top) == len(bottom), "EDTA: rows differ in length" assert without_gaps(top) == s, "EDTA: top row does not spell s" assert without_gaps(bottom) == t, "EDTA: bottom row does not spell t" assert mismatches == distance, "EDTA: alignment costs " + str(mismatches) + ", not " + str(distance) } ``` ## CTEA — Counting Optimal Alignments [Problem statement](https://rosalind.info/problems/ctea/) ```biolang # Rosalind: CTEA — Counting Optimal Alignments # https://rosalind.info/problems/ctea/ # # Given: Two protein strings s and t. # Return: The number of optimal alignments between them, modulo 134,217,727. let s = "PLEASANTLY" let t = "MEANLY" let modulus = 134217727 fn set_row(table, index, row) { range(0, len(table)) |> map(|i| if i == index then row else table[i]) } let rows = len(s) + 1 let cols = len(t) + 1 # First the edit distances, exactly as in EDTA. let cost = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 then j else i)) let i = 1 while i < rows { let si = substr(s, i - 1, 1) let row = [i] let j = 1 while j < cols { let step = if si == substr(t, j - 1, 1) then 0 else 1 row = push(row, min([cost[i - 1][j - 1] + step, cost[i - 1][j] + 1, row[j - 1] + 1])) j = j + 1 } cost = set_row(cost, i, row) i = i + 1 } # Then count the paths that achieve those distances. A cell's count is the sum # of the counts of every predecessor that reaches it at optimal cost, so ties # multiply out rather than being collapsed. let ways = range(0, rows) |> map(|i| range(0, cols) |> map(|j| if i == 0 or j == 0 then 1 else 0)) let a = 1 while a < rows { let sa = substr(s, a - 1, 1) let row = [1] let b = 1 while b < cols { let step = if sa == substr(t, b - 1, 1) then 0 else 1 let best = cost[a][b] # Bound separately: a continued expression cannot start a line with `+`. let via_diagonal = if cost[a - 1][b - 1] + step == best then ways[a - 1][b - 1] else 0 let via_up = if cost[a - 1][b] + 1 == best then ways[a - 1][b] else 0 let via_left = if cost[a][b - 1] + 1 == best then row[b - 1] else 0 let total = via_diagonal + via_up + via_left row = push(row, total % modulus) b = b + 1 } ways = set_row(ways, a, row) a = a + 1 } let distance = cost[rows - 1][cols - 1] let result = ways[rows - 1][cols - 1] println("Edit distance: " + str(distance)) println("Result: " + str(result)) println("Expected: 4 optimal alignments at distance 5") fn test_ctea_optimal_alignment_count() { assert distance == 5, "CTEA: distance " + str(distance) assert result == 4, "CTEA: got " + str(result) # There is always at least one optimal alignment. assert result >= 1, "CTEA: counted none" } ``` ## FOUN — The Founder Effect and Genetic Drift [Problem statement](https://rosalind.info/problems/foun/) ```biolang # Rosalind: FOUN — The Founder Effect and Genetic Drift # https://rosalind.info/problems/foun/ # # Given: Positive integers N and m, and an array A of m allele counts. # Return: An N x m matrix whose (i, j) entry is the common logarithm of the # probability that the recessive allele is lost entirely after i+1 generations, # starting from A[j] copies in a population of N diploid individuals. let n = 4 let counts = [0, 1, 2] let alleles = 2 * n fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) } fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) } # One Wright-Fisher generation: from i copies, the next count is # Binomial(2N, i / 2N) — every allele is drawn afresh. fn step(distribution) { range(0, alleles + 1) |> map(|j| { range(0, alleles + 1) |> map(|i| { let p = float(i) / float(alleles) distribution[i] * float(choose(alleles, j)) * pow(p, float(j)) * pow(1.0 - p, float(alleles - j)) }) |> sum() }) } # Each row is one further generation; each column one starting count. let rows = counts |> map(|start| { let distribution = range(0, alleles + 1) |> map(|i| if i == start then 1.0 else 0.0) range(0, n) |> map(|_| { distribution = step(distribution) log10(distribution[0]) }) }) println("Result (rows are starting counts, columns are generations):") range(0, len(counts)) |> each(|j| { println(" from " + str(counts[j]) + ": " + (rows[j] |> map(|v| str(round(v, 6))) |> join(" "))) }) println("") println("Expected, generation 1: 0.0 -0.463935 -0.999510") println(" from 1 copy: (7/8)^8 = 0.343609, log10 = -0.463936") println(" from 2 copies: (6/8)^8 = 0.100113, log10 = -0.999510") fn test_foun_founder_effect() { # An allele already absent stays absent, so the probability is 1 and its # logarithm 0 in every generation. # Bound first: `x |> count_if(f) == n` reads the comparison as count_if's # second argument rather than as a test of its result. let certain = rows[0] |> count_if(|v| v == 0.0) assert certain == n, "FOUN: the zero-copy column is not certain" assert round(rows[1][0], 6) == -0.463936, "FOUN: one copy, one generation = " + str(rows[1][0]) assert round(rows[2][0], 6) == -0.99951, "FOUN: two copies, one generation = " + str(rows[2][0]) # Loss can only become more likely with time, so each row is non-decreasing. let i = 1 while i < n { assert rows[1][i] >= rows[1][i - 1], "FOUN: loss became less likely at generation " + str(i) i = i + 1 } } ``` ## GCON — Global Alignment with Constant Gap Penalty [Problem statement](https://rosalind.info/problems/gcon/) ```biolang # Rosalind: GCON — Global Alignment with Constant Gap Penalty # https://rosalind.info/problems/gcon/ # # Given: Two protein strings. # Return: The maximum global alignment score using BLOSUM62 and a constant gap # penalty of 5 — a run of gaps costs 5 however long it is. let s = "PLEASANTLY" let t = "MEANLY" let gap_open = -5.0 let gap_extend = 0.0 let blosum = score_matrix("blosum62") let residues = blosum.row_names let width = blosum.ncol fn residue_index(names, ch) { (range(0, len(names)) |> filter(|i| names[i] == ch))[0] } fn substitution(mat, names, cols, a, b) { mat.data[residue_index(names, a) * cols + residue_index(names, b)] } # The affine recurrence with a zero extension cost: opening a gap is charged # once and continuing it is free, which is exactly "constant". fn constant_gap_score(a, b, mat, names, cols, open_penalty, extend_penalty) { let n = len(b) let very_low = -1000000.0 let prev_m = concat([0.0], range(1, n + 1) |> map(|_| very_low)) let prev_x = concat([very_low], range(1, n + 1) |> map(|_| very_low)) let prev_y = concat([very_low], range(1, n + 1) |> map(|_| open_penalty)) let i = 1 while i <= len(a) { let ai = substr(a, i - 1, 1) let row_m = [very_low] let row_x = [open_penalty] let row_y = [very_low] let j = 1 while j <= n { let sub = substitution(mat, names, cols, ai, substr(b, j - 1, 1)) let best_prev = max([prev_m[j - 1], prev_x[j - 1], prev_y[j - 1]]) row_m = push(row_m, best_prev + sub) row_x = push(row_x, max([prev_m[j] + open_penalty, prev_x[j] + extend_penalty])) row_y = push(row_y, max([row_m[j - 1] + open_penalty, row_y[j - 1] + extend_penalty])) j = j + 1 } prev_m = row_m prev_x = row_x prev_y = row_y i = i + 1 } max([prev_m[n], prev_x[n], prev_y[n]]) } let result = constant_gap_score(s, t, blosum, residues, width, gap_open, gap_extend) println("Result: " + str(int(result))) println("Expected: 13 — the same pair scores 8 under GLOB's per-residue gap of 5,") println(" so charging the four-residue gap once instead of four times") println(" recovers 3 x 5 = 15 minus the 10 already counted.") fn test_gcon_constant_gap_alignment() { assert int(result) == 13, "GCON: got " + str(result) # A constant gap penalty can never score worse than the linear one, which # charges the same opening plus more for every extra residue. assert int(result) >= 8, "GCON: scored below the linear-gap result" } ``` ## PRSM — Matching a Spectrum to a Protein [Problem statement](https://rosalind.info/problems/prsm/) ```biolang # Rosalind: PRSM — Matching a Spectrum to a Protein # https://rosalind.info/problems/prsm/ # # Given: A collection of protein strings and a multiset R of masses. # Return: The largest multiplicity of any protein's complete spectrum convolved # with R, and a protein achieving it. let proteins = ["GSDMQS", "VWICN", "IASMQS", "PVSMGAD"] let spectrum = [445.17838, 115.02694, 186.07931, 314.13789, 317.1198, 215.09061] let monoisotopic = { A: 71.03711, C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841, G: 57.02146, H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406, M: 131.04049, N: 114.04293, P: 97.05276, Q: 128.05858, R: 156.10111, S: 87.03203, T: 101.04768, V: 99.06841, W: 186.07931, Y: 163.06333 } fn mass_of(peptide) { range(0, len(peptide)) |> reduce(|acc, i| acc + monoisotopic[substr(peptide, i, 1)], 0.0) } # The complete spectrum is every prefix and every suffix mass — the fragments a # spectrometer would see if the peptide broke at each bond in turn. fn complete_spectrum(peptide) { let prefixes = range(1, len(peptide)) |> map(|k| mass_of(substr(peptide, 0, k))) let suffixes = range(1, len(peptide)) |> map(|k| mass_of(substr(peptide, len(peptide) - k, k))) concat(prefixes, suffixes) } # Convolving the two multisets, the winning shift is the parent mass offset; # its multiplicity is how many fragments line up. fn best_multiplicity(peptide) { let own = complete_spectrum(peptide) let shifts = spectrum |> flat_map(|r| own |> map(|m| round(r - m, 5))) let tallies = shifts |> unique() |> map(|d| shifts |> count_if(|x| x == d)) max(tallies) } let scored = proteins |> map(|p| { protein: p, score: best_multiplicity(p) }) let best = scored |> sort_by(|e| e.score) |> reverse() |> first() println("Scores:") scored |> each(|e| println(" " + e.protein + ": " + str(e.score))) println("Result: " + str(best.score)) println(" " + best.protein) println("Expected: 3 and IASMQS") println("") println("Note: GSDMQS also reaches 3, so the maximum is shared. Rosalind accepts") println("any protein achieving it; which one is reported here depends on the") println("order of the input rather than on anything meaningful.") let maximal = scored |> filter(|e| e.score == best.score) |> map(|e| e.protein) fn test_prsm_spectrum_to_protein() { assert best.score == 3, "PRSM: multiplicity " + str(best.score) # The tie is real, so require the winner to be one of the maximal proteins # rather than one particular string. assert maximal |> contains(best.protein), "PRSM: winner is not maximal" assert maximal |> contains("IASMQS"), "PRSM: IASMQS is not among the maximal proteins" # A complete spectrum holds two fragments per bond. assert len(complete_spectrum("IASMQS")) == 2 * (len("IASMQS") - 1), "PRSM: spectrum size is wrong" } ``` ## PDPL — Creating a Restriction Map [Problem statement](https://rosalind.info/problems/pdpl/) ```biolang # Rosalind: PDPL — Creating a Restriction Map # https://rosalind.info/problems/pdpl/ # # Given: A multiset L containing every pairwise difference of a set of n # points on a line, together with 0. # Return: A set of points whose pairwise differences are exactly L. let differences = [0, 2, 2, 3, 3, 4, 5, 6, 7, 8, 10] # The largest difference spans the whole map, so the first and last points are # fixed at once. Everything else is placed by trying each remaining distance as # the next point and checking it does not claim a difference L cannot supply — # the classic turnpike reconstruction. let width = max(differences) let candidates = differences |> filter(|d| d > 0 and d < width) |> unique() |> sort() fn pairwise(points) { range(0, len(points)) |> flat_map(|i| { range(i + 1, len(points)) |> map(|j| points[j] - points[i]) }) |> sort() } let target = differences |> filter(|d| d > 0) |> sort() # n points produce n(n-1)/2 differences, which fixes how many to place. let n = len(target) let point_count = 2 while point_count * (point_count - 1) / 2 < n { point_count = point_count + 1 } let interior_needed = point_count - 2 # Every combination of interior points, smallest first, until one reproduces L. fn combinations(items, size) { if size == 0 { [[]] } else { range(0, len(items)) |> flat_map(|i| { combinations(drop(items, i + 1), size - 1) |> map(|rest| concat([items[i]], rest)) }) } } let solution = (combinations(candidates, interior_needed) |> map(|interior| concat(concat([0], interior), [width])) |> filter(|points| pairwise(points) == target))[0] println("Result: " + (solution |> map(|p| str(p)) |> join(" "))) println("Expected: 0 2 4 7 10") println("Check: its pairwise differences are " + (pairwise(solution) |> map(|d| str(d)) |> join(" "))) fn test_pdpl_restriction_map() { assert len(solution) == point_count, "PDPL: placed " + str(len(solution)) + " points" # The answer is only correct if its differences reproduce the input exactly, # which is a stronger check than matching one printed string. assert pairwise(solution) == target, "PDPL: differences do not match L" assert solution[0] == 0, "PDPL: does not start at 0" assert solution[len(solution) - 1] == width, "PDPL: does not end at the widest span" } ``` ## CSET — Fixing an Inconsistent Character Set [Problem statement](https://rosalind.info/problems/cset/) ```biolang # Rosalind: CSET — Fixing an Inconsistent Character Set # https://rosalind.info/problems/cset/ # # Given: A collection of characters over the same taxa, at most one of which is # inconsistent with the others. # Return: The remaining characters, which form a consistent set. let characters = ["100001", "000110", "111000", "100111"] # Two splits are compatible when one of the four ways of intersecting their # halves is empty — the four-gamete condition. If all four appear, no tree can # carry both characters at once. fn compatible(a, b) { let seen = range(0, len(a)) |> map(|i| substr(a, i, 1) ++ substr(b, i, 1)) |> unique() len(seen) < 4 } fn all_compatible(rows) { let bad = range(0, len(rows)) |> flat_map(|i| { range(i + 1, len(rows)) |> filter(|j| !compatible(rows[i], rows[j])) }) len(bad) == 0 } # Drop each character in turn and keep the first set that becomes consistent. let kept = (range(0, len(characters)) |> map(|i| concat(take(characters, i), drop(characters, i + 1))) |> filter(|rows| all_compatible(rows)))[0] let removed = characters |> filter(|c| !(kept |> contains(c))) println("Result:") kept |> each(|c| println(" " + c)) println("Removed: " + (removed |> join(" "))) println("Expected: one character removed, leaving a consistent set") println(" (100001 and 111000 are the incompatible pair here)") fn test_cset_consistent_character_set() { assert len(kept) == len(characters) - 1, "CSET: kept " + str(len(kept)) + " characters" # The point of the answer is consistency, so check it rather than a string. assert all_compatible(kept), "CSET: the kept set is still inconsistent" # The original set really was inconsistent, otherwise the problem is empty. assert !all_compatible(characters), "CSET: the input was already consistent" } ``` ## LREP — Finding the Longest Multiple Repeat [Problem statement](https://rosalind.info/problems/lrep/) ```biolang # Rosalind: LREP — Finding the Longest Multiple Repeat # https://rosalind.info/problems/lrep/ # # Given: A DNA string s with $ appended, a positive integer k, and the edges of # the suffix tree of s. # Return: The longest substring of s occurring at least k times. let s = "CATACATAC$" let k = 2 # The suffix tree comes with the problem, so nothing has to build one. Each edge # gives its parent, its child, where the label starts in s (1-based) and how long # it is. let tree_edges = [ { parent: "node1", child: "node2", start: 1, length: 1 }, { parent: "node1", child: "node7", start: 2, length: 1 }, { parent: "node1", child: "node14", start: 3, length: 3 }, { parent: "node1", child: "node17", start: 10, length: 1 }, { parent: "node2", child: "node3", start: 2, length: 4 }, { parent: "node2", child: "node6", start: 10, length: 1 }, { parent: "node3", child: "node4", start: 6, length: 5 }, { parent: "node3", child: "node5", start: 10, length: 1 }, { parent: "node7", child: "node8", start: 3, length: 3 }, { parent: "node7", child: "node11", start: 5, length: 1 }, { parent: "node8", child: "node9", start: 6, length: 5 }, { parent: "node8", child: "node10", start: 10, length: 1 }, { parent: "node11", child: "node12", start: 6, length: 5 }, { parent: "node11", child: "node13", start: 10, length: 1 }, { parent: "node14", child: "node15", start: 6, length: 5 }, { parent: "node14", child: "node16", start: 10, length: 1 }, ] # A substring's number of occurrences is the number of leaves below the node it # ends at, because every leaf is one suffix that starts with it. So the answer is # the deepest node with at least k leaves under it, and the string is the path # from the root spelled out along the way. let children = {} let has_children = {} for edge in tree_edges { if contains(keys(children), edge.parent) { children[edge.parent] = push(children[edge.parent], edge) } else { children[edge.parent] = [edge] } has_children[edge.child] = true } # Walk down from the root, carrying the label spelled so far. A node with no # children is a leaf and counts as one occurrence; anything else is the sum of # its children. Written as an explicit stack because the recursion would have to # return two things at once. let best = "" let stack = [{ node: "node1", label: "" }] let leaves = {} # First pass: how many leaves hang below each node, deepest first. let order = [] while len(stack) > 0 { let top = stack[len(stack) - 1] stack = slice(stack, 0, len(stack) - 1) order = push(order, top) if contains(keys(children), top.node) { for edge in children[top.node] { let label = top.label + substr(s, edge.start - 1, edge.length) stack = push(stack, { node: edge.child, label: label }) } } } # `order` holds parents before children, so counting it backwards means every # child is already counted when its parent is reached. let i = len(order) - 1 while i >= 0 { let entry = order[i] if contains(keys(children), entry.node) { let total = 0 for edge in children[entry.node] { total = total + leaves[edge.child] } leaves[entry.node] = total # The root spells nothing, and a repeat has to be a real substring. if total >= k and len(entry.label) > len(best) { best = entry.label } } else { leaves[entry.node] = 1 } i = i - 1 } println("Result: " + best) println("Expected: CATAC") fn test_lrep_longest_multiple_repeat() { assert best == "CATAC", "LREP: got " + best # It has to occur at least k times, and be a substring of s. assert len(find_motif(dna(substr(s, 0, len(s) - 1)), dna(best))) >= k, "LREP: " + best + " does not occur " + str(k) + " times" } ``` ## MREP — Identifying Maximal Repeats [Problem statement](https://rosalind.info/problems/mrep/) ```biolang # Rosalind: MREP — Identifying Maximal Repeats # https://rosalind.info/problems/mrep/ # # Given: A DNA string s of length at most 1 kbp. # Return: Every maximal repeat of s of length at least 20. # # A repeat is maximal when it occurs at least twice and no pair of its # occurrences can be extended by one symbol in either direction and still agree. let s = "TAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTATTATATAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTAT" let minimum = 20 # Sort the suffixes, then read the repeats off the LCP array: lcp[i] is how much # suffix i shares with the one before it, so any value >= 20 is a candidate # repeat, and the two suffixes it relates are two of its occurrences. # # The sentinel keeps a suffix that is a prefix of another from being extended # past the end of the string. let text = s + "$" let sa = suffix_array(text) let lcp = lcp_array(text) # Right-maximal: a candidate is right-maximal when it is not simply the start of # a longer shared prefix, which the LCP array says directly — take each candidate # at its own length rather than a prefix of it. # # Left-maximal: the characters immediately before the two occurrences must # differ. If every occurrence is preceded by the same symbol, the repeat extends # leftwards and is not maximal. fn left_character(position) { if position == 0 then "^" else substr(text, position - 1, 1) } let found = [] for i in range(1, len(sa)) { let length = lcp[i] if length >= minimum { let candidate = substr(text, sa[i], length) # Already recorded from an earlier pair. if contains(found, candidate) == false { # Collect every occurrence, then ask whether they all share a left # neighbour. Positions come from the original string, not the # sentinel-terminated one. let occurrences = find_motif(dna(s), dna(candidate)) let lefts = occurrences |> map(|p| left_character(p)) |> unique() if len(occurrences) >= 2 and len(lefts) > 1 { found = push(found, candidate) } } } } # The longest first, which is how Rosalind's sample reads. let result = found |> sort_by(|a, b| len(b) - len(a)) println("Result:") for repeat_seq in result { println(" " + repeat_seq) } println("Expected:") println(" TAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTAT") println(" ATGGGTCCAGAGTTTTGTAATTT") fn test_mrep_maximal_repeats() { assert len(result) == 2, "MREP: expected 2 maximal repeats, got " + str(len(result)) assert contains(result, "TAGAGATAGAATGGGTCCAGAGTTTTGTAATTTCCATGGGTCCAGAGTTTTGTAATTTAT"), "MREP: the longer repeat is missing from " + str(result) assert contains(result, "ATGGGTCCAGAGTTTTGTAATTT"), "MREP: the shorter repeat is missing from " + str(result) } ``` ## LAFF — Local Alignment with Affine Gap Penalty [Problem statement](https://rosalind.info/problems/laff/) ```biolang # Rosalind: LAFF — Local Alignment with Affine Gap Penalty # https://rosalind.info/problems/laff/ # # Given: Two protein strings s and t. # Return: The maximum local alignment score, and the substrings of s and t that # achieve it. BLOSUM62, gap opening 11, gap extension 1. let s = "PLEASANTLY" let t = "MEANLY" # Local mode forgives both ends and lets the alignment stop early; the affine # gap is the gap_open/gap_extend pair. A gap of length L costs 11 + (L-1), so # opening is -10 on top of the -1 every symbol pays. let result = align(s, t, "local", 0, 0, -1, -10, "blosum62") # Local alignment reports the aligned substrings; strip the gaps to recover the # pieces of s and t themselves. fn without_gaps(aligned) { range(0, len(aligned)) |> map(|i| substr(aligned, i, 1)) |> filter(|c| c != "-") |> join("") } let sub_s = without_gaps(result.aligned_a) let u = without_gaps(result.aligned_b) println("Result: " + str(result.score)) println(" " + sub_s) println(" " + u) println("Expected: 12 / LEAS / MEAN") fn test_laff_local_affine_alignment() { assert result.score == 12, "LAFF: got " + str(result.score) assert contains(s, sub_s), "LAFF: " + sub_s + " is not a substring of s" assert contains(t, u), "LAFF: " + u + " is not a substring of t" } ``` ## SMGB — Semiglobal Alignment [Problem statement](https://rosalind.info/problems/smgb/) ```biolang # Rosalind: SMGB — Semiglobal Alignment # https://rosalind.info/problems/smgb/ # # Given: Two DNA strings s and t. # Return: The maximum semiglobal alignment score, and an alignment achieving it. # Match +1, substitution -1, linear gap -1. let s = dna"CAGCACTTGGATTCTCGG" let t = dna"CAGCGTGG" # Semiglobal forgives gaps at both ends of both sequences, so neither overhang # is charged for — which is what lets a short sequence sit inside a long one # without paying for the flanks. Global alignment of these two scores far worse. let result = align(s, t, "semiglobal", 1, -1, -1, 0) let global_score = align(s, t, "global", 1, -1, -1, 0).score println("Result: " + str(result.score)) println("Expected: 4") println("Alignment:") println(" " + result.aligned_a) println(" " + result.aligned_b) println("(global alignment of the same pair scores " + str(global_score) + ")") fn test_smgb_semiglobal_alignment() { assert result.score == 4, "SMGB: got " + str(result.score) # Forgiving the end gaps can only help, never hurt. assert result.score >= global_score, "SMGB: semiglobal scored below global" } ``` ## MGAP — Maximizing the Gap Symbols of an Optimal Alignment [Problem statement](https://rosalind.info/problems/mgap/) ```biolang # Rosalind: MGAP — Maximizing the Gap Symbols of an Optimal Alignment # https://rosalind.info/problems/mgap/ # # Given: Two DNA strings s and t. # Return: The largest number of gap symbols that can appear in any maximum-score # alignment of s and t, for any scoring with match > 0 and both penalties < 0. let s = dna"AACGTA" let t = dna"ACACCTA" # The scoring is left open on purpose, and that is the whole problem: since a # match is worth something and everything else costs, a maximum-score alignment # must match as many symbols as it possibly can. That is the longest common # subsequence. Every symbol outside it is opposite a gap — len(s) - lcs of them # in one row and len(t) - lcs in the other. let common = lcs(str(s), str(t)) let gaps = len(str(s)) + len(str(t)) - 2 * len(common) println("LCS: " + common) println("Result: " + str(gaps)) println("Expected: 3") fn test_mgap_maximum_gap_symbols() { assert gaps == 3, "MGAP: got " + str(gaps) # The subsequence has to be a real one of both strings. assert is_subsequence(common, str(s)), "MGAP: LCS is not a subsequence of s" assert is_subsequence(common, str(t)), "MGAP: LCS is not a subsequence of t" } ``` ## CSTR — Creating a Character Table from Genetic Strings [Problem statement](https://rosalind.info/problems/cstr/) ```biolang # Rosalind: CSTR — Creating a Character Table from Genetic Strings # https://rosalind.info/problems/cstr/ # # Given: A collection of characterizable DNA strings of equal length. # Return: A character table whose rows are the nontrivial characters, one per # position where the strings disagree. let strings = [ "ATGCTACC", "CGTTTACC", "ATTCGACC", "AGTCTCCC", "CGTCTATC", ] let width = len(strings[0]) # A position splits the collection into two groups. The split is nontrivial only # when both groups have at least two members: a position where one string # differs from all the others separates nothing, and neither does a position # where every string agrees. let characters = [] for column in range(0, width) { let symbols = strings |> map(|s| substr(s, column, 1)) let distinct = unique(symbols) if len(distinct) == 2 { let first_count = symbols |> filter(|c| c == distinct[0]) |> len() let second_count = len(symbols) - first_count if first_count >= 2 and second_count >= 2 { # Which state gets 1 is arbitrary, so take the first symbol seen. characters = push(characters, symbols |> map(|c| if c == distinct[0] then "1" else "0") |> join("")) } } } println("Result:") for character in characters { println(" " + character) } println("Expected:") println(" 10110") println(" 10100") fn test_cstr_character_table() { assert len(characters) == 2, "CSTR: expected 2 nontrivial characters, got " + str(len(characters)) # The 0/1 assignment is arbitrary, so a row and its complement are the same # character. Compare each row against the expected one either way round. fn flipped(row) { range(0, len(row)) |> map(|i| if substr(row, i, 1) == "1" then "0" else "1") |> join("") } assert characters[0] == "10110" or flipped(characters[0]) == "10110", "CSTR: first row " + characters[0] assert characters[1] == "10100" or flipped(characters[1]) == "10100", "CSTR: second row " + characters[1] } ``` ## SUFF — Encoding Suffix Trees [Problem statement](https://rosalind.info/problems/suff/) ```biolang # Rosalind: SUFF — Encoding Suffix Trees # https://rosalind.info/problems/suff/ # # Given: A DNA string s ending in $. # Return: The substrings labelling the edges of the suffix tree of s. let text = "ATAAATG$" # Built from the suffix array and its LCP array rather than by inserting every # suffix into a trie and collapsing it. Both give the same tree; the difference # is cost. A trie of all suffixes has O(n^2) nodes before anything is collapsed, # whereas the two arrays are linear in the text and the tree can be read off them # in a single left-to-right pass. # # That is why aligners store the arrays and never materialise the tree: for a # human genome the tree does not fit, and the arrays do. let sa = suffix_array(text) let lcp = lcp_array(text) # A leaf hangs from the deeper of what its suffix shares with either neighbour, # so its edge is whatever is left of the suffix below that depth. let leaf_labels = range(0, len(sa)) |> map(|i| { let before = if i == 0 then 0 else lcp[i] let after = if i + 1 < len(sa) then lcp[i + 1] else 0 let depth = max([before, after]) substr(text, sa[i] + depth, len(text) - sa[i] - depth) }) # Internal nodes come from the LCP array alone: a run of suffixes sharing a # prefix of length L hangs off one node at depth L. A stack tracks which are # still open as the array is scanned. let internal_labels = [] let stack = [{ depth: 0, start: 0 }] for i in range(1, len(sa)) { let shared = lcp[i] let last_start = sa[i] while stack[len(stack) - 1].depth > shared { let closing = stack[len(stack) - 1] stack = slice(stack, 0, len(stack) - 1) let parent_depth = max([stack[len(stack) - 1].depth, shared]) internal_labels = push(internal_labels, substr(text, closing.start + parent_depth, closing.depth - parent_depth)) last_start = closing.start } if stack[len(stack) - 1].depth < shared { stack = push(stack, { depth: shared, start: last_start }) } } while len(stack) > 1 { let closing = stack[len(stack) - 1] stack = slice(stack, 0, len(stack) - 1) let parent_depth = stack[len(stack) - 1].depth internal_labels = push(internal_labels, substr(text, closing.start + parent_depth, closing.depth - parent_depth)) } let labels = leaf_labels + internal_labels println("Result: " + join(sort(labels), " ")) println("Expected (any order): AAATG$ G$ T ATG$ TG$ A A AAATG$ G$ T G$ $") fn test_suff_encoding_suffix_trees() { let expected = ["AAATG$", "G$", "T", "ATG$", "TG$", "A", "A", "AAATG$", "G$", "T", "G$", "$"] assert sort(labels) == sort(expected), "SUFF: got " + join(sort(labels), " ") # One leaf per suffix, each reaching the sentinel. assert len(leaf_labels) == len(text), "SUFF: one leaf per suffix" for label in leaf_labels { assert substr(label, len(label) - 1, 1) == "$", "SUFF: a leaf edge must reach the end" } # Every label is a real substring, and the whole tree spells every suffix. for label in labels { assert contains(text, label), "SUFF: " + label + " is not in the text" } # Concatenating the edges costs less than storing the suffixes outright — # which is the saving a suffix tree exists for. let edge_characters = labels |> map(|l| len(l)) |> sum() let all_suffixes = range(0, len(text)) |> map(|i| len(text) - i) |> sum() assert edge_characters < all_suffixes, "SUFF: the tree should be smaller than the suffixes it encodes" } ``` ## SGRA — Using the Spectrum Graph to Infer Peptides [Problem statement](https://rosalind.info/problems/sgra/) ```biolang # Rosalind: SGRA — Using the Spectrum Graph to Infer Peptides # https://rosalind.info/problems/sgra/ # # Given: A list of positive real numbers (a spectrum). # Return: The longest protein string matching the spectrum graph. let spectrum = [ 3524.8542, 3623.5245, 3710.9335, 3841.974, 3929.00603, 3970.0326, 4026.05879, 4057.0646, 4083.08025, ] # Monoisotopic masses, not the integer table the cyclopeptide problems use. # Real instrument readings carry decimals, so a difference is matched within a # tolerance rather than exactly — 1e-4 here, which is far tighter than the gap # between any two residue masses and far looser than float noise. let monoisotopic = { A: 71.03711, C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841, G: 57.02146, H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406, M: 131.04049, N: 114.04293, P: 97.05276, Q: 128.05858, R: 156.10111, S: 87.03203, T: 101.04768, V: 99.06841, W: 186.07931, Y: 163.06333 } let tolerance = 0.0001 # Where L and I collide exactly, and K and Q nearly do, take the first — a # spectrum cannot separate them, so naming either is equally correct. let residues = ["G", "A", "S", "P", "V", "T", "C", "L", "N", "D", "Q", "K", "E", "M", "H", "F", "R", "Y", "W"] fn residue_between(lighter, heavier, table, letters, slack) { let gap = heavier - lighter let hits = letters |> filter(|letter| abs(table[letter] - gap) < slack) if len(hits) > 0 then hits[0] else "" } # The graph is a DAG — masses only increase — so the longest path is found by # scanning in sorted order and never revisiting a node. Searching the paths # themselves would be exponential; here each node is settled once, from the best # way of reaching it. let ordered = sort(spectrum) let best = ordered |> map(|_| "") # Where each best path began, so the answer's mass can be checked against the # span it was actually read across. The path skips nodes, so this cannot be # recovered by counting backwards from the end. let origin = range(0, len(ordered)) for j in range(0, len(ordered)) { for i in range(0, j) { let letter = residue_between(ordered[i], ordered[j], monoisotopic, residues, tolerance) if letter != "" and len(best[i]) + 1 > len(best[j]) { best[j] = best[i] + letter origin[j] = origin[i] } } } let finish = argmax(best |> map(|s| len(s))) let answer = best[finish] println("Result: " + answer) println("Expected: WMSPG") fn test_sgra_spectrum_graph() { assert answer == "WMSPG", "SGRA: got " + answer # Every residue of the answer must be a real gap between two spectrum values. assert len(answer) == 5, "SGRA: expected a 5-residue peptide" # The peptide's own mass is the span it was read across. let peptide_mass = chars(answer) |> map(|c| monoisotopic[c]) |> sum() let span = ordered[finish] - ordered[origin[finish]] assert abs(peptide_mass - span) < len(answer) * tolerance, "SGRA: the peptide weighs " + str(round(peptide_mass, 4)) + " but spans " + str(round(span, 4)) # Nothing longer exists in the graph. for candidate in best { assert len(candidate) <= len(answer), "SGRA: a longer path was missed" } } ``` ## FULL — Inferring Peptide from Full Spectrum [Problem statement](https://rosalind.info/problems/full/) ```biolang # Rosalind: FULL — Inferring Peptide from Full Spectrum # https://rosalind.info/problems/full/ # # Given: 2n+3 positive reals — a parent mass, then the b-ions and y-ions of a # peptide of length n, in no particular order. # Return: A protein string of length n consistent with them. let parent_mass = 1988.21104821 let ions = [ 610.391039105, 738.485999105, 766.492149105, 863.544909105, 867.528589105, 992.587499105, 995.623549105, 1120.6824591, 1124.6661391, 1221.7188991, 1249.7250491, 1377.8200091, ] let monoisotopic = { A: 71.03711, C: 103.00919, D: 115.02694, E: 129.04259, F: 147.06841, G: 57.02146, H: 137.05891, I: 113.08406, K: 128.09496, L: 113.08406, M: 131.04049, N: 114.04293, P: 97.05276, Q: 128.05858, R: 156.10111, S: 87.03203, T: 101.04768, V: 99.06841, W: 186.07931, Y: 163.06333 } let residues = ["G", "A", "S", "P", "V", "T", "C", "L", "N", "D", "Q", "K", "E", "M", "H", "F", "R", "Y", "W"] let tolerance = 0.001 # A peptide fragments from both ends at once. A b-ion is a prefix, a y-ion the # matching suffix, so the two always sum to the parent mass — and the list mixes # them with no labels saying which is which. # # The trick is that it does not matter. Walk upwards from the lightest ion; the # next ion differing from it by a residue mass is the next prefix, whichever kind # it nominally is. Taking an ion also spends its complement, since a prefix and # its suffix are one fragmentation event and cannot both extend the chain. let ordered = sort(ions) let used = ordered |> map(|_| false) fn residue_between(lighter, heavier, table, letters, slack) { let gap = heavier - lighter let hits = letters |> filter(|letter| abs(table[letter] - gap) < slack) if len(hits) > 0 then hits[0] else "" } fn index_of_mass(masses, wanted, slack) { let hits = range(0, len(masses)) |> filter(|i| abs(masses[i] - wanted) < slack) if len(hits) > 0 then hits[0] else 0 - 1 } fn spend(marks, masses, at, total, slack) { marks[at] = true let partner = index_of_mass(masses, total - masses[at], slack) if partner >= 0 { marks[partner] = true } marks } let peptide = "" let current = 0 used = spend(used, ordered, current, parent_mass, tolerance) let searching = true while searching { let onward = range(current + 1, len(ordered)) |> filter(|j| used[j] == false and residue_between(ordered[current], ordered[j], monoisotopic, residues, tolerance) != "") if len(onward) == 0 { searching = false } else { let next_index = onward[0] peptide = peptide + residue_between(ordered[current], ordered[next_index], monoisotopic, residues, tolerance) used = spend(used, ordered, next_index, parent_mass, tolerance) current = next_index } } println("Result: " + peptide) println("Expected: KEKEP") fn test_full_inferring_peptide_from_full_spectrum() { assert peptide == "KEKEP", "FULL: got " + peptide # n is fixed by the input's size, and the answer has to match it. let n = (len(ions) + 1 - 3) / 2 assert len(peptide) == n, "FULL: expected " + str(n) + " residues" # Every ion pairs with another summing to the parent mass — the property that # makes b-ions and y-ions indistinguishable here, and the reason taking one # spends the other. for mass in ordered { let partner = index_of_mass(ordered, parent_mass - mass, tolerance) assert partner >= 0, "FULL: " + str(round(mass, 4)) + " has no complement" } assert len(ordered) % 2 == 0, "FULL: the ions must pair up" } ``` ## ALPH — Alignment-Based Phylogeny [Problem statement](https://rosalind.info/problems/alph/) ```biolang # Rosalind: ALPH — Alignment-Based Phylogeny # https://rosalind.info/problems/alph/ # # Given: A rooted binary tree in Newick format and a multiple alignment of its # leaves. # Return: The minimum total Hamming distance over the tree's edges, and internal # node labels achieving it. # (((ostrich,cat)rat,(duck,fly)mouse)dog,(elephant,pikachu)hamster)robot; let children = { "rat": ["ostrich", "cat"], "mouse": ["duck", "fly"], "dog": ["rat", "mouse"], "hamster": ["elephant", "pikachu"], "robot": ["dog", "hamster"], } let root = "robot" let leaves = { "ostrich": "AC", "cat": "CA", "duck": "T-", "fly": "GC", "elephant": "-T", "pikachu": "AA", } # The gap is a fifth symbol here, not a missing value. That is the whole reason # this problem is separate from the ordinary parsimony ones: in an alignment a # gap is evidence — an insertion or deletion actually happened — and treating it # as unknown would make every column containing one free. let symbols = ["A", "C", "G", "T", "-"] let width = len(leaves["ostrich"]) fn is_leaf(node, leaf_map) { contains(keys(leaf_map), node) } # Sankoff: the cheapest cost of a subtree given its root's symbol is the sum over # children of their cheapest cost plus one if the symbol has to change. Columns # are independent, so each is solved separately and the scores added. fn score_column(node, position, leaf_map, child_map, alphabet) { if is_leaf(node, leaf_map) { let here = substr(leaf_map[node], position, 1) return { costs: alphabet |> map(|s| if s == here then 0 else 1000000), picks: {}, solved: [], kids: [], } } let kids = child_map[node] let solved = kids |> map(|kid| score_column(kid, position, leaf_map, child_map, alphabet)) let picks = {} let costs = range(0, len(alphabet)) |> map(|mine| { range(0, len(kids)) |> map(|k| { let options = range(0, len(alphabet)) |> map(|theirs| solved[k].costs[theirs] + (if theirs == mine then 0 else 1)) let chosen = argmin(options) picks[str(mine) + "," + str(k)] = chosen options[chosen] }) |> sum() }) { costs: costs, picks: picks, solved: solved, kids: kids } } fn assign(node, tree, chosen_index, labels, alphabet) { labels[node] = labels[node] + alphabet[chosen_index] for k in range(0, len(tree.kids)) { if len(tree.solved[k].kids) > 0 { labels = assign(tree.kids[k], tree.solved[k], tree.picks[str(chosen_index) + "," + str(k)], labels, alphabet) } } labels } let labels = {} for node in keys(children) { labels[node] = "" } let total = 0 for position in range(0, width) { let solved = score_column(root, position, leaves, children, symbols) let best = argmin(solved.costs) total = total + solved.costs[best] labels = assign(root, solved, best, labels, symbols) } let named = {} for node in keys(leaves) { named[node] = leaves[node] } for node in keys(children) { named[node] = labels[node] } fn hamming(a, b) { range(0, len(a)) |> count_if(|i| substr(a, i, 1) != substr(b, i, 1)) } println("Result: " + str(total)) for node in sort(keys(children)) { println(" " + node + ": " + named[node]) } println("Expected: 8 (rat AC, mouse TC, dog AC, hamster AT, robot AC)") println("This labelling differs and costs the same — several are optimal, and the") println("assertion recounts the changes rather than trusting the printed one.") fn test_alph_alignment_based_phylogeny() { assert total == 8, "ALPH: scored " + str(total) # The score has to be the changes the labelling actually shows — a number # without a matching labelling is the usual bug here. let shown = sort(keys(children)) |> flat_map(|parent| children[parent] |> map(|kid| hamming(named[parent], named[kid]))) |> sum() assert shown == total, "ALPH: the labelling shows " + str(shown) + " changes but the score is " + str(total) # Every internal label is the right length and drawn from the alphabet. for node in keys(children) { assert len(named[node]) == width, "ALPH: " + node + " has the wrong length" assert (chars(named[node]) |> count_if(|c| contains(symbols, c) == false)) == 0, "ALPH: " + named[node] + " uses a symbol outside the alphabet" } # Treating the gap as free would score lower, which is exactly what must not # happen — it is a real evolutionary event. assert (chars(leaves["duck"]) |> count_if(|c| c == "-")) == 1, "ALPH: duck carries a gap" } ``` ## GASM — Genome Assembly Using Reads [Problem statement](https://rosalind.info/problems/gasm/) ```biolang # Rosalind: GASM — Genome Assembly Using Reads # https://rosalind.info/problems/gasm/ # # Given: Error-free reads of equal length, whose de Bruijn graph is two directed # cycles. # Return: A cyclic superstring of minimal length containing every read or its # reverse complement. let reads = ["AATCT", "TGTAA", "GATTA", "ACAGA"] # A read gives no clue which strand it came from, so the assembly graph has to # hold every read *and* its reverse complement — which is why the de Bruijn graph # here always falls into exactly two cycles, one the mirror of the other. Either # spells the answer; they are the same circular genome read from opposite # strands. # # k is not given. The right one is the largest that still leaves every node with # somewhere to go, so it is searched for downwards from the read length: too # large and the graph breaks into fragments, too small and distinct repeats # collapse into one node. let both_strands = reads |> flat_map(|r| [r, str(reverse_complement(dna(r)))]) # Every k-length window of every read, not the reads themselves. The reads are # longer than the k that works, and using them whole leaves the graph with more # nodes than edges — which is what a broken assembly looks like. fn windows_of(patterns, k) { patterns |> flat_map(|p| range(0, len(p) - k + 1) |> map(|i| substr(p, i, k))) |> unique() } fn cycle_through(k, patterns) { let edges_of = windows_of(patterns, k) let graph = {} for pattern in edges_of { let prefix = substr(pattern, 0, k - 1) let suffix = substr(pattern, 1, k - 1) if contains(keys(graph), prefix) { graph[prefix] = push(graph[prefix], suffix) } else { graph[prefix] = [suffix] } } # Exactly one way in and one way out of every node, or this k does not give a # clean pair of cycles. let nodes = unique(keys(graph) + (keys(graph) |> flat_map(|n| graph[n]))) let single_exit = nodes |> count_if(|n| contains(keys(graph), n) and len(graph[n]) == 1) if single_exit != len(nodes) { return "" } let start = sort(keys(graph))[0] let walk = [start] let at = graph[start][0] while at != start { walk = push(walk, at) at = graph[at][0] } # Half the nodes belong to the mirror cycle, so a correct walk covers exactly # half of them. if len(walk) * 2 != len(nodes) { return "" } # Circular, so only the leading character of each node is kept — the trailing # k-1 are the leading ones come round again. walk |> map(|node| substr(node, 0, 1)) |> join("") } let answer = "" let k = len(reads[0]) while answer == "" and k > 2 { answer = cycle_through(k, both_strands) k = k - 1 } println("Result: " + answer) println("Expected: GATTACA (any rotation, of either strand, is accepted)") fn test_gasm_genome_assembly() { assert len(answer) == 7, "GASM: expected a 7-character superstring, got " + str(len(answer)) # Every read, or its reverse complement, must appear when the string is read # around the circle. let wrapped = answer + substr(answer, 0, len(reads[0]) - 1) for read in reads { let flipped = str(reverse_complement(dna(read))) assert contains(wrapped, read) or contains(wrapped, flipped), "GASM: neither " + read + " nor " + flipped + " occurs" } # GATTACA is one rotation of one strand; the answer is correct up to both. let rotations = range(0, len(answer)) |> map(|i| substr(wrapped, i, len(answer))) let mirror = str(reverse_complement(dna(answer))) let mirror_wrapped = mirror + substr(mirror, 0, len(answer)) let mirror_rotations = range(0, len(answer)) |> map(|i| substr(mirror_wrapped, i, len(answer))) assert contains(rotations, "GATTACA") or contains(mirror_rotations, "GATTACA"), "GASM: " + answer + " is not GATTACA up to rotation and strand" } ``` ## MULT — Multiple Alignment [Problem statement](https://rosalind.info/problems/mult/) ```biolang # Rosalind: MULT — Multiple Alignment # https://rosalind.info/problems/mult/ # # Given: Four DNA strings of length at most 10. # Return: A multiple alignment of maximum score, where every mismatched pair in # a column costs 1 and matched symbols — including two gaps — cost nothing. let sequences = ["ATATCCG", "TCCG", "ATGTACTG", "ATGTCTG"] # Two sequences need a table, three a cube, four a hypercube — and the number of # moves out of each cell is 2^k - 1, every non-empty choice of which sequences # advance. Fifteen here, against three for a pairwise alignment. # # That is the whole lesson of this problem: exact multiple alignment costs # O(n^k), so it stops being possible at a handful of short sequences. Every tool # that aligns hundreds of sequences is approximating, usually by aligning pairs # and merging. let moves = range(1, 16) |> map(|mask| range(0, 4) |> map(|s| (mask / pow(2, s)) % 2 |> floor())) let lengths = sequences |> map(|s| len(s)) fn key(at) { at |> map(|v| str(v)) |> join(",") } # A column costs one for every disagreeing pair. Two gaps agree; a gap against a # base does not. fn column_cost(symbols) { range(0, len(symbols)) |> flat_map(|i| range(i + 1, len(symbols)) |> filter(|j| symbols[i] != symbols[j])) |> len() } let best = {} let came_from = {} best[key([0, 0, 0, 0])] = 0 for i in range(0, lengths[0] + 1) { for j in range(0, lengths[1] + 1) { for k in range(0, lengths[2] + 1) { for l in range(0, lengths[3] + 1) { let here = [i, j, k, l] if i + j + k + l > 0 { let at = key(here) best[at] = 0 - 1000000 for move in moves { let previous = range(0, 4) |> map(|s| here[s] - move[s]) if (previous |> count_if(|v| v < 0)) == 0 { let symbols = range(0, 4) |> map(|s| if move[s] == 1 then substr(sequences[s], previous[s], 1) else "-") let candidate = best[key(previous)] - column_cost(symbols) if candidate > best[at] { best[at] = candidate came_from[at] = move } } } } } } } } let score = best[key(lengths)] let rows = ["", "", "", ""] let at = lengths while (at |> sum()) > 0 { let move = came_from[key(at)] for s in range(0, 4) { if move[s] == 1 { rows[s] = substr(sequences[s], at[s] - 1, 1) + rows[s] at[s] = at[s] - 1 } else { rows[s] = "-" + rows[s] } } } println("Result: " + str(score)) for line in rows { println(" " + line) } println("Expected: -18") println(" ATAT-CCG / -T---CCG / ATGTACTG / ATGT-CTG") fn test_mult_multiple_alignment() { assert score == 0 - 18, "MULT: scored " + str(score) # Structural checks, since several alignments reach the optimum. let widths = rows |> map(|r| len(r)) |> unique() assert len(widths) == 1, "MULT: every row must be the same length" for s in range(0, 4) { assert replace(rows[s], "-", "") == sequences[s], "MULT: row " + str(s) + " is not its original sequence" } # The alignment shown must actually score what was claimed. let recounted = 0 - (range(0, widths[0]) |> map(|c| column_cost(rows |> map(|r| substr(r, c, 1)))) |> sum()) assert recounted == score, "MULT: the alignment shown scores " + str(recounted) + ", not " + str(score) # No column is entirely gaps — that move does not exist, and one would be # free under this scoring. let empty_columns = range(0, widths[0]) |> count_if(|c| (rows |> count_if(|r| substr(r, c, 1) == "-")) == 4) assert empty_columns == 0, "MULT: an all-gap column would be free and meaningless" } ``` ## GREP — Genome Assembly with Perfect Coverage and Repeats [Problem statement](https://rosalind.info/problems/grep/) ```biolang # Rosalind: GREP — Genome Assembly with Perfect Coverage and Repeats # https://rosalind.info/problems/grep/ # # Given: (k+1)-mers from one strand of a circular chromosome. # Return: Every circular string assembled by a complete cycle in the de Bruijn # graph, each beginning with the first (k+1)-mer given. let patterns = [ "CAG", "AGT", "GTT", "TTT", "TTG", "TGG", "GGC", "GCG", "CGT", "GTT", "TTC", "TCA", "CAA", "AAT", "ATT", "TTC", "TCA", ] # PCOV assumes perfect coverage and returns one answer. This is the honest # version: repeats make the Eulerian cycle non-unique, and every cycle is a # genome consistent with the reads. Six here — so the reads simply do not # determine the chromosome, and reporting one would be picking arbitrarily. # # That ambiguity is why repeats are the hard part of real assembly, and why # BA3J's paired reads exist: extra distance information cuts the alternatives # down. let k = len(patterns[0]) - 1 let edges_list = patterns |> map(|p| { from_node: substr(p, 0, k), to_node: substr(p, 1, k), }) let start_node = edges_list[0].from_node # Depth-first over the edge multiset. Hierholzer finds one cycle in linear time; # finding them all is a search, and only feasible because the graph is small. let complete = [] let stack = [{ at: edges_list[0].to_node, used: range(0, len(edges_list)) |> map(|i| i == 0), walk: [start_node, edges_list[0].to_node], }] while len(stack) > 0 { let state = stack[len(stack) - 1] stack = slice(stack, 0, len(stack) - 1) let remaining = state.used |> count_if(|u| u == false) if remaining == 0 { # A cycle only counts if it closed. if state.at == start_node { complete = push(complete, state.walk) } } else { for i in range(0, len(edges_list)) { if state.used[i] == false and edges_list[i].from_node == state.at { let marked = state.used marked[i] = true stack = push(stack, { at: edges_list[i].to_node, used: marked, walk: push(state.walk, edges_list[i].to_node), }) } } } } # A circular string keeps one character per edge — the trailing k are the leading # k come round again. let assembled = complete |> map(|walk| range(0, len(walk) - 1) |> map(|i| substr(walk[i], 0, 1)) |> join("")) |> unique() |> sort() println("Result:") for line in assembled { println(" " + line) } println("Expected: CAGTTCAATTTGGCGTT CAGTTCAATTGGCGTTT CAGTTTCAATTGGCGTT") println(" CAGTTTGGCGTTCAATT CAGTTGGCGTTCAATTT CAGTTGGCGTTTCAATT") fn test_grep_assembly_with_repeats() { let expected = ["CAGTTCAATTTGGCGTT", "CAGTTCAATTGGCGTTT", "CAGTTTCAATTGGCGTT", "CAGTTTGGCGTTCAATT", "CAGTTGGCGTTCAATTT", "CAGTTGGCGTTTCAATT"] assert assembled == sort(expected), "GREP: got " + join(assembled, " ") # Each answer is as long as there are reads, and starts with the first one. for answer in assembled { assert len(answer) == len(patterns), "GREP: " + answer + " should be " + str(len(patterns)) + " long" assert substr(answer, 0, k + 1) == patterns[0], "GREP: " + answer + " must begin with " + patterns[0] # And read around the circle, its composition is exactly the input. let wrapped = answer + substr(answer, 0, k) let composition = range(0, len(answer)) |> map(|i| substr(wrapped, i, k + 1)) assert sort(composition) == sort(patterns), "GREP: " + answer + " does not have the given composition" } # The point of the problem: more than one genome fits. assert len(assembled) > 1, "GREP: repeats should leave the assembly ambiguous" } ``` ## MPRT — Finding a Protein Motif [Problem statement](https://rosalind.info/problems/mprt/) ```biolang # @skip (needs a network connection — remove this line to run it) # Rosalind: MPRT — Finding a Protein Motif # https://rosalind.info/problems/mprt/ # # Given: UniProt access IDs. # Return: For each protein containing the N-glycosylation motif, its ID and the # 1-based positions where the motif occurs. # # This example fetches from UniProt, so it runs in the advisory job rather than # the hermetic gate every other Stronghold problem passes. let ids = ["A2Z669", "B5ZC00", "P07204_TRBM_HUMAN", "P20840_SAG1_YEAST"] # N{P}[ST]{P} — asparagine, then anything but proline, then serine or threonine, # then anything but proline. Not a fixed string, which is the point: a motif is a # pattern with alternatives and exclusions, and searching for it is not # substring matching. fn has_motif_at(protein, i) { substr(protein, i, 1) == "N" and substr(protein, i + 1, 1) != "P" and (substr(protein, i + 2, 1) == "S" or substr(protein, i + 2, 1) == "T") and substr(protein, i + 3, 1) != "P" } fn motif_positions(protein) { range(0, len(protein) - 3) |> filter(|i| has_motif_at(protein, i)) |> map(|i| i + 1) } # UniProt is keyed by the accession alone; the trailing name in an ID like # P07204_TRBM_HUMAN is Rosalind's own annotation. fn accession_of(id) { split(id, "_")[0] } fn sequence_of(id) { let fasta = str(uniprot_fasta(accession_of(id))) lines(fasta) |> filter(|line| substr(line, 0, 1) != ">") |> join("") } let found = ids |> map(|id| { id: id, positions: motif_positions(sequence_of(id)) }) |> filter(|entry| len(entry.positions) > 0) println("Result:") for entry in found { println(" " + entry.id) println(" " + (entry.positions |> map(|p| str(p)) |> join(" "))) } println("Expected: B5ZC00 85 118 142 306 395") println(" P07204_TRBM_HUMAN 47 115 116 382 409") println(" P20840_SAG1_YEAST 79 109 135 248 306 348 364 402 485 501 614") fn test_mprt_finding_a_protein_motif() { let expected = { "B5ZC00": [85, 118, 142, 306, 395], "P07204_TRBM_HUMAN": [47, 115, 116, 382, 409], "P20840_SAG1_YEAST": [79, 109, 135, 248, 306, 348, 364, 402, 485, 501, 614], } assert len(found) == 3, "MPRT: expected 3 proteins with the motif, got " + str(len(found)) for entry in found { assert contains(keys(expected), entry.id), "MPRT: unexpected protein " + entry.id assert entry.positions == expected[entry.id], "MPRT: " + entry.id + " gave " + str(entry.positions) } # A2Z669 has no motif, and leaving it out of the answer is part of the task. assert (found |> count_if(|e| e.id == "A2Z669")) == 0, "MPRT: A2Z669 has no motif" # The matcher itself, checked without the network: the exclusions are what # make this a motif rather than a substring search. assert motif_positions("NASA") == [1], "MPRT: NAS_ matches" assert motif_positions("NPSA") == [], "MPRT: proline in the second position blocks it" assert motif_positions("NASP") == [], "MPRT: proline in the fourth blocks it" assert motif_positions("NAGA") == [], "MPRT: the third must be S or T" assert motif_positions("NATA") == [1], "MPRT: threonine works as well as serine" } ``` ## OSYM — Isolating Symbols in Alignments [Problem statement](https://rosalind.info/problems/osym/) ```biolang # Rosalind: OSYM — Isolating Symbols in Alignments # https://rosalind.info/problems/osym/ # # Given: Two DNA strings. # Return: The maximum global alignment score, and the sum over all (j,k) of the # best score of any alignment that pairs s[j] with t[k]. # Matches score +1, mismatches and gaps -1. let first_string = "ATAGATA" let second_string = "ACAGGTA" let gap = 0 - 1 fn substitution(a, b) { if a == b then 1 else 0 - 1 } # The question is not what the best alignment is, but how good the best # alignment *through each particular pairing* would be. Recomputing an alignment # for every one of the 49 pairs would be quadratic work repeated quadratically. # # Instead: the best alignment through (j,k) is the best way of reaching it plus # the best way of leaving it. One table filled forwards and one filled backwards # answer every pair at once — the same trick BA5K uses to find a middle edge. fn forward_table(rows, columns, indel) { let table = [range(0, len(columns) + 1) |> map(|j| j * indel)] for i in range(1, len(rows) + 1) { let line = [i * indel] for j in range(1, len(columns) + 1) { let diagonal = table[i - 1][j - 1] + substitution(substr(rows, i - 1, 1), substr(columns, j - 1, 1)) let up = table[i - 1][j] + indel let left = line[j - 1] + indel line = push(line, max([diagonal, up, left])) } table = push(table, line) } table } let forward = forward_table(first_string, second_string, gap) # The backward table is the forward one on both strings reversed, so only one # direction has to be written. let backward_reversed = forward_table(reverse(first_string), reverse(second_string), gap) let n = len(first_string) let m = len(second_string) fn backward_at(table, i, j, rows, columns) { # Score of aligning s[i..) with t[j..), read out of the reversed table. table[rows - i][columns - j] } let best_score = forward[n][m] let total = range(1, n + 1) |> flat_map(|j| range(1, m + 1) |> map(|k| { let paired = substitution(substr(first_string, j - 1, 1), substr(second_string, k - 1, 1)) forward[j - 1][k - 1] + paired + backward_at(backward_reversed, j, k, n, m) })) |> sum() println("Result: " + str(best_score)) println(" " + str(total)) println("Expected: 3") println(" -139") fn test_osym_isolating_symbols() { assert best_score == 3, "OSYM: the best alignment scores " + str(best_score) assert total == 0 - 139, "OSYM: the sum is " + str(total) # No pairing can beat the unconstrained optimum, and at least one must reach # it — the best alignment pairs something. let every = range(1, n + 1) |> flat_map(|j| range(1, m + 1) |> map(|k| { let paired = substitution(substr(first_string, j - 1, 1), substr(second_string, k - 1, 1)) forward[j - 1][k - 1] + paired + backward_at(backward_reversed, j, k, n, m) })) assert max(every) == best_score, "OSYM: the best constrained score should equal the unconstrained one" for value in every { assert value <= best_score, "OSYM: a constrained alignment cannot beat the optimum" } assert len(every) == n * m, "OSYM: one entry per pair of positions" } ``` ## ITWV — Finding Disjoint Motifs in a Gene [Problem statement](https://rosalind.info/problems/itwv/) ```biolang # Rosalind: ITWV — Finding Disjoint Motifs in a Gene # https://rosalind.info/problems/itwv/ # # Given: A DNA string s and a collection of patterns. # Return: The matrix M where M[j][k] = 1 if patterns j and k can be interwoven # into some substring of s. let text = "GACCACGGTT" let patterns = ["ACAG", "GT", "CCG"] # Two motifs are interwoven when a stretch of the sequence spells both at once, # each as a subsequence, using every character for at most one of them. That is # what makes this different from finding each motif separately: they compete for # the same characters, so the question is whether a shared stretch can satisfy # both — which matters when asking whether two binding sites can overlap. # # The state is (how much of t is placed, how much of u is placed), and the # characters of s are consumed in order. Each character extends whichever motif # it matches, so a cell is reachable from at most two others. fn can_interweave(source, start, left, right) { let rows = len(left) let columns = len(right) # reached[a][b] — the first a of `left` and first b of `right` are placed. let reached = range(0, rows + 1) |> map(|_| range(0, columns + 1) |> map(|_| false)) reached[0] = reached[0] let seed = reached[0] seed[0] = true reached[0] = seed let at = start let done = false while at < len(source) and done == false { let symbol = substr(source, at, 1) let next = range(0, rows + 1) |> map(|a| range(0, columns + 1) |> map(|_| false)) for a in range(0, rows + 1) { let line = next[a] for b in range(0, columns + 1) { if reached[a][b] { # The character can be skipped only by ending the window, so # every reachable state must consume it or the window stops # here. Carrying it forward unconsumed is what would let the # two motifs be found in unrelated parts of the sequence. if a < rows and substr(left, a, 1) == symbol { let advanced = next[a + 1] advanced[b] = true next[a + 1] = advanced } if b < columns and substr(right, b, 1) == symbol { line[b + 1] = true } } } next[a] = line } reached = next if reached[rows][columns] { done = true } at = at + 1 } done } fn interwoven_anywhere(source, left, right) { (range(0, len(source)) |> count_if(|start| can_interweave(source, start, left, right))) > 0 } let pairings = range(0, len(patterns)) |> map(|j| range(0, len(patterns)) |> map(|k| if interwoven_anywhere(text, patterns[j], patterns[k]) then 1 else 0)) println("Result:") for line in pairings { println(" " + (line |> map(|v| str(v)) |> join(" "))) } println("Expected: 0 0 1 / 0 1 0 / 1 0 0") fn test_itwv_disjoint_motifs() { let shown = pairings |> map(|line| line |> map(|v| str(v)) |> join(" ")) |> join(" / ") assert shown == "0 0 1 / 0 1 0 / 1 0 0", "ITWV: got " + shown # The relation is symmetric — interweaving t with u is the same question as # interweaving u with t. for j in range(0, len(patterns)) { for k in range(0, len(patterns)) { assert pairings[j][k] == pairings[k][j], "ITWV: the pairings must be symmetric" } } # GT with itself works: GACCACGGTT contains GGTT, which spells GT twice using # disjoint characters. assert pairings[1][1] == 1, "ITWV: GT interweaves with itself" # ACAG with itself does not — there are not enough characters to spell it # twice disjointly anywhere. assert pairings[0][0] == 0, "ITWV: ACAG cannot be interwoven with itself" } ``` ## CNTQ — Counting Quartets [Problem statement](https://rosalind.info/problems/cntq/) ```biolang # Rosalind: CNTQ — Counting Quartets # https://rosalind.info/problems/cntq/ # # Given: n and an unrooted binary tree on n taxa in Newick format. # Return: The number of quartets consistent with the tree, modulo 1,000,000. let taxa = ["lobster", "cat", "dog", "caterpillar", "elephant", "mouse"] # (lobster,(cat,dog),(caterpillar,(elephant,mouse))); # Each internal edge splits the taxa in two; these are the non-trivial sides. let clades = [ ["cat", "dog"], ["caterpillar", "elephant", "mouse"], ["elephant", "mouse"], ] # The answer is simply C(n,4), and the reason is worth more than the number: in a # *fully resolved* binary tree, any four taxa are separated by some internal edge # into two pairs, so every four-taxon set contributes exactly one consistent # quartet. Nothing about the tree's shape enters into it. # # That stops being true the moment the tree has an unresolved node, which is why # QRT and CHBP — where splits are partial — are real problems and this one is # not. let n = len(taxa) let answer = choose(n, 4) % 1000000 println("Result: " + str(answer)) println("Expected: 15") fn test_cntq_counting_quartets() { assert answer == 15, "CNTQ: got " + str(answer) # Demonstrated rather than asserted: enumerate all four-taxon subsets and # confirm each really is separated into two pairs by some edge. fn subsets_of_four(items) { range(0, len(items)) |> flat_map(|a| range(a + 1, len(items)) |> flat_map(|b| range(b + 1, len(items)) |> flat_map(|c| range(c + 1, len(items)) |> map(|d| [items[a], items[b], items[c], items[d]])))) } let quartets = subsets_of_four(taxa) assert len(quartets) == choose(n, 4), "CNTQ: C(6,4) is 15 subsets" let consistent = quartets |> count_if(|four| (clades |> count_if(|clade| (four |> count_if(|t| contains(clade, t))) == 2)) > 0) assert consistent == len(quartets), "CNTQ: " + str(consistent) + " of " + str(len(quartets)) + " subsets are separated" assert consistent == answer, "CNTQ: the count and the formula must agree" } ``` ## QRT — Quartets [Problem statement](https://rosalind.info/problems/qrt/) ```biolang # Rosalind: QRT — Quartets # https://rosalind.info/problems/qrt/ # # Given: A partial character table. # Return: Every quartet inferable from the splits its characters describe. let taxa = ["cat", "dog", "elephant", "ostrich", "mouse", "rabbit", "robot"] let characters = [ "01xxx00", "x11xx00", "111x00x", ] # A partial character is one that could not be scored for every taxon — an 'x' # means "not known", not "third state". So a character constrains only the taxa # it actually saw, and the quartets it supports are pairs drawn from its two # scored groups. # # The whole point is that missing data does not make a character useless. It # still separates the taxa it did score, and those separations are what a tree # gets built from. fn group_of(character, symbol, names) { range(0, len(names)) |> filter(|i| substr(character, i, 1) == symbol) |> map(|i| names[i]) } fn pairs_of(items) { range(0, len(items)) |> flat_map(|a| range(a + 1, len(items)) |> map(|b| [items[a], items[b]])) } # A quartet has no orientation — {a,b}|{c,d} is the same as {c,d}|{a,b} — so each # is recorded in a canonical form and duplicates dropped. Two characters often # support the same quartet, and counting it twice would overstate the evidence. fn canonical(left, right) { let one = join(sort(left), ",") let two = join(sort(right), ",") if one < two then one + "|" + two else two + "|" + one } let seen = {} let quartets = [] for character in characters { let zeros = group_of(character, "0", taxa) let ones = group_of(character, "1", taxa) for left in pairs_of(zeros) { for right in pairs_of(ones) { let key = canonical(left, right) if contains(keys(seen), key) == false { seen[key] = true quartets = push(quartets, { left: left, right: right }) } } } } let written = quartets |> map(|q| "{" + join(q.left, ", ") + "} {" + join(q.right, ", ") + "}") |> sort() println("Result:") for line in written { println(" " + line) } println("Expected (any order): {elephant, dog} {rabbit, robot} / {cat, dog} {mouse, rabbit}") println(" {mouse, rabbit} {cat, elephant} / {dog, elephant} {mouse, rabbit}") fn test_qrt_quartets() { assert len(quartets) == 4, "QRT: expected 4 quartets, got " + str(len(quartets)) # Compared as unordered pairs of unordered pairs, since neither side nor the # order within a side carries meaning. let expected = [ canonical(["elephant", "dog"], ["rabbit", "robot"]), canonical(["cat", "dog"], ["mouse", "rabbit"]), canonical(["mouse", "rabbit"], ["cat", "elephant"]), canonical(["dog", "elephant"], ["mouse", "rabbit"]), ] let got = quartets |> map(|q| canonical(q.left, q.right)) assert sort(got) == sort(expected), "QRT: got " + join(sort(got), " ") # Every quartet's four taxa are distinct and really were scored by some # character — an 'x' can never appear in one. for q in quartets { let four = q.left + q.right assert len(unique(four)) == 4, "QRT: a quartet must name four different taxa" let supporting = characters |> count_if(|c| (q.left |> count_if(|t| substr(c, (range(0, len(taxa)) |> filter(|i| taxa[i] == t))[0], 1) == "0")) == 2 and (q.right |> count_if(|t| substr(c, (range(0, len(taxa)) |> filter(|i| taxa[i] == t))[0], 1) == "1")) == 2) assert supporting > 0, "QRT: no character supports " + canonical(q.left, q.right) } } ``` ## CHBP — Character-Based Phylogeny [Problem statement](https://rosalind.info/problems/chbp/) ```biolang # Rosalind: CHBP — Character-Based Phylogeny # https://rosalind.info/problems/chbp/ # # Given: Species names and a consistent character table. # Return: An unrooted binary tree in Newick format modelling the table. let taxa = ["cat", "dog", "elephant", "mouse", "rabbit", "rat"] let characters = ["011101", "001101", "001100"] # The inverse of CTBL: there, a tree produced splits; here the splits have to # produce the tree. It works because a consistent table's splits are *laminar* — # any two are nested or disjoint, never crossing — and a laminar family is # exactly a tree. # # An unrooted tree has no first node, so a reference taxon is picked and every # split is taken from the side away from it. That makes the clades nest, and the # reference becomes the outermost branch. let reference = taxa[0] fn side_away_from(character, names, anchor) { let anchor_index = (range(0, len(names)) |> filter(|i| names[i] == anchor))[0] let anchor_state = substr(character, anchor_index, 1) range(0, len(names)) |> filter(|i| substr(character, i, 1) != anchor_state) |> map(|i| names[i]) } let clades = characters |> map(|c| side_away_from(c, taxa, reference)) fn is_subset(small, big) { (small |> count_if(|x| contains(big, x) == false)) == 0 } fn build(members, groups) { # The groups strictly inside this one, and of those the maximal — anything # contained in another is handled a level deeper. let inside = groups |> filter(|g| is_subset(g, members) and len(g) < len(members)) let maximal = inside |> filter(|g| (inside |> count_if(|other| len(other) > len(g) and is_subset(g, other))) == 0) let covered = maximal |> flat_map(|g| g) let loose = members |> filter(|m| contains(covered, m) == false) let parts = (maximal |> map(|g| build(g, groups))) + loose if len(parts) == 1 then parts[0] else "(" + join(parts, ",") + ")" } let rest = taxa |> filter(|t| t != reference) let inner = build(rest, clades) # The reference sits alongside the rest at the unrooted centre, so its branch is # spliced in rather than wrapped around. let newick = "(" + reference + "," + substr(inner, 1, len(inner) - 2) + ");" println("Result: " + newick) println("Expected: (dog,(cat,rabbit),(rat,(elephant,mouse)));") println("Both describe one unrooted tree — they differ only in which branch is") println("written first, which an unrooted tree does not fix.") fn test_chbp_character_based_phylogeny() { # The real requirement is that the tree induces exactly the table's splits. # Newick is not canonical for an unrooted tree, so comparing strings would be # comparing a formatting choice. fn splits_of(tree_text, names) { # Every parenthesised group is a clade; take each as a split. let found = [] let stack = [] for i in range(0, len(tree_text)) { let symbol = substr(tree_text, i, 1) if symbol == "(" { stack = push(stack, i) } if symbol == ")" { let opened = stack[len(stack) - 1] stack = slice(stack, 0, len(stack) - 1) let inner_text = substr(tree_text, opened + 1, i - opened - 1) let members = names |> filter(|t| contains(inner_text, t)) if len(members) > 1 and len(members) < len(names) { found = push(found, join(sort(members), ",")) } } } unique(found) } let mine = splits_of(newick, taxa) let published = splits_of("(dog,(cat,rabbit),(rat,(elephant,mouse)));", taxa) # A split and its complement are the same split, so compare canonically. fn canonical_splits(raw, names) { raw |> map(|s| { let members = split(s, ",") let other = names |> filter(|t| contains(members, t) == false) let one = join(sort(members), ",") let two = join(sort(other), ",") if one < two then one else two }) |> unique() |> sort() } assert canonical_splits(mine, taxa) == canonical_splits(published, taxa), "CHBP: splits " + join(canonical_splits(mine, taxa), " | ") + " against " + join(canonical_splits(published, taxa), " | ") # And those splits are exactly the ones the characters describe. let from_characters = canonical_splits( clades |> map(|c| join(sort(c), ",")), taxa) assert canonical_splits(mine, taxa) == from_characters, "CHBP: the tree must induce the table's splits and no others" } ``` ## EUBT — Enumerating Unrooted Binary Trees [Problem statement](https://rosalind.info/problems/eubt/) ```biolang # Rosalind: EUBT — Enumerating Unrooted Binary Trees # https://rosalind.info/problems/eubt/ # # Given: n taxa. # Return: Every unrooted binary tree on those taxa, in Newick format. let taxa = ["dog", "cat", "mouse", "elephant"] # Built by insertion: start with the only tree on three taxa — a star — then add # each remaining taxon by splitting one existing edge in two and hanging it off # the new node. Every unrooted binary tree arises exactly once this way, which is # why the count is the double factorial (2n-5)!! and why it explodes: fifteen # trees at n=5, over two million at n=10. # # That growth is the reason nobody enumerates trees to find the best one. It is # also why the parsimony and distance methods in this pack exist. let edges_start = [ { a: "center", b: taxa[0] }, { a: "center", b: taxa[1] }, { a: "center", b: taxa[2] }, ] let trees = [edges_start] for step in range(3, len(taxa)) { let leaf = taxa[step] let grown = [] for tree in trees { for i in range(0, len(tree)) { let fresh = "node" + str(step) + "_" + str(i) let split_edge = tree[i] let others = range(0, len(tree)) |> filter(|j| j != i) |> map(|j| tree[j]) grown = push(grown, others + [ { a: split_edge.a, b: fresh }, { a: fresh, b: split_edge.b }, { a: fresh, b: leaf }, ]) } } trees = grown } fn neighbours_of(tree, node) { (tree |> filter(|e| e.a == node) |> map(|e| e.b)) + (tree |> filter(|e| e.b == node) |> map(|e| e.a)) } fn is_taxon(node, names) { contains(names, node) } fn walk(tree, node, parent, names) { if is_taxon(node, names) { return node } let children = neighbours_of(tree, node) |> filter(|c| c != parent) "(" + (children |> map(|c| walk(tree, c, node, names)) |> join(",")) + ")" } # Written from the first taxon, which an unrooted tree does not privilege — it # is just somewhere to start reading. fn to_newick(tree, names) { let anchor = names[0] let hub = neighbours_of(tree, anchor)[0] let branches = neighbours_of(tree, hub) |> filter(|c| c != anchor) "(" + anchor + "," + (branches |> map(|c| walk(tree, c, hub, names)) |> join(",")) + ");" } let written = trees |> map(|t| to_newick(t, taxa)) println("Result:") for line in written { println(" " + line) } println("Expected: three trees — (mouse,cat)|(elephant,dog), (elephant,mouse)|(cat,dog),") println(" (elephant,cat)|(mouse,dog), written from a different anchor.") fn test_eubt_enumerating_unrooted_binary_trees() { # (2n-5)!! trees: 1 x 3 for n = 4. assert len(trees) == 3, "EUBT: expected 3 trees, got " + str(len(trees)) assert len(unique(written)) == 3, "EUBT: the trees must be distinct" # A tree on n taxa has 2n-3 edges and n-2 internal nodes. for tree in trees { assert len(tree) == 2 * len(taxa) - 3, "EUBT: expected " + str(2 * len(taxa) - 3) + " edges, got " + str(len(tree)) for taxon in taxa { assert len(neighbours_of(tree, taxon)) == 1, "EUBT: " + taxon + " must be a leaf" } } # The three topologies are exactly the three ways of pairing four taxa. That # is the whole content of the answer, and it does not depend on how the # Newick is rooted. fn pairing_of(tree, names) { # The one internal edge separates the taxa into two pairs. Which side the # edge happens to be stored from is arbitrary, so the smaller of the two # names the topology. let internal = tree |> filter(|e| is_taxon(e.a, names) == false and is_taxon(e.b, names) == false) let side = neighbours_of(tree, internal[0].a) |> filter(|x| is_taxon(x, names)) let other = names |> filter(|t| contains(side, t) == false) let one = join(sort(side), ",") let two = join(sort(other), ",") if one < two then one else two } let pairings = trees |> map(|t| pairing_of(t, taxa)) |> sort() assert len(unique(pairings)) == 3, "EUBT: each tree must pair the taxa differently" # With four taxa the topology is fixed by which one cat is paired with, and # all three possibilities appear exactly once. assert pairings == sort(["cat,dog", "cat,mouse", "cat,elephant"]), "EUBT: got pairings " + join(pairings, " | ") } ``` ## QRTD — Quartet Distance [Problem statement](https://rosalind.info/problems/qrtd/) ```biolang # Rosalind: QRTD — Quartet Distance # https://rosalind.info/problems/qrtd/ # # Given: n taxa and two unrooted binary trees. # Return: The quartet distance between them. let taxa = ["A", "B", "C", "D", "E"] # (A,C,((B,D),E)); and (C,(B,D),(A,E)); # Each internal edge is a split; these are the non-trivial sides. let first_clades = [["B", "D"], ["B", "D", "E"]] let second_clades = [["B", "D"], ["A", "E"]] # Two trees can share every taxon and still disagree about almost everything, so # comparing them needs a measure. Quartets give one: each set of four taxa is # separated into two pairs by exactly one edge, and two trees either agree about # that pairing or they do not. # # It beats counting shared splits because it degrades gracefully — moving one # taxon changes a handful of quartets, where it can destroy every split at once. fn pairing_in(four, clades) { # The clade cutting these four 2-and-2 names the quartet. let cutting = clades |> filter(|c| (four |> count_if(|t| contains(c, t))) == 2) if len(cutting) == 0 { return "" } let side = four |> filter(|t| contains(cutting[0], t)) let other = four |> filter(|t| contains(cutting[0], t) == false) let one = join(sort(side), ",") let two = join(sort(other), ",") if one < two then one + "|" + two else two + "|" + one } fn subsets_of_four(items) { range(0, len(items)) |> flat_map(|a| range(a + 1, len(items)) |> flat_map(|b| range(b + 1, len(items)) |> flat_map(|c| range(c + 1, len(items)) |> map(|d| [items[a], items[b], items[c], items[d]])))) } let quartets = subsets_of_four(taxa) let shared = quartets |> count_if(|four| { let one = pairing_in(four, first_clades) let two = pairing_in(four, second_clades) one != "" and one == two }) # Both trees are fully resolved, so each induces exactly one quartet per subset. let first_count = quartets |> count_if(|four| pairing_in(four, first_clades) != "") let second_count = quartets |> count_if(|four| pairing_in(four, second_clades) != "") let distance = first_count + second_count - 2 * shared println("Result: " + str(distance)) println("Expected: 4") println("(" + str(len(quartets)) + " quartets each, " + str(shared) + " agreeing)") fn test_qrtd_quartet_distance() { assert distance == 4, "QRTD: got " + str(distance) # A resolved tree resolves every four-taxon subset, which CNTQ established. assert first_count == choose(len(taxa), 4), "QRTD: the first tree resolves all 5" assert second_count == choose(len(taxa), 4), "QRTD: and so does the second" assert shared == 3, "QRTD: 3 of the 5 quartets agree" # A tree against itself is distance zero — the identity any distance must # satisfy, and worth checking rather than assuming. let self_shared = quartets |> count_if(|four| pairing_in(four, first_clades) != "") assert first_count + first_count - 2 * self_shared == 0, "QRTD: a tree must be distance 0 from itself" # The two disagreeing quartets are the ones involving both B,D and A,E. let disagreeing = quartets |> filter(|four| pairing_in(four, first_clades) != pairing_in(four, second_clades)) assert len(disagreeing) == 2, "QRTD: exactly two quartets differ" } ``` ## RNAS — Wobble Bonding and RNA Secondary Structures [Problem statement](https://rosalind.info/problems/rnas/) ```biolang # 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" } ``` ## KSIM — Finding All Similar Motifs [Problem statement](https://rosalind.info/problems/ksim/) ```biolang # Rosalind: KSIM — Finding All Similar Motifs # https://rosalind.info/problems/ksim/ # # Given: k, a motif s, and a genome t. # Return: Every substring of t within edit distance k of s, as (position, length). let k = 2 let motif = "ACGTAG" let sequence = "ACGGATCGGCATCGT" # Approximate matching, which is what real motif search always is: a binding site # is a tendency rather than a fixed string, and an exact search finds a fraction # of the real sites. # # A caveat worth stating rather than hiding: this checks every (start, length) # pair directly, which is O(n^2) edit-distance computations. It is clear and it # is correct, and at 15 bases it is instant — but the problem permits a 50 kbp # genome, where it would not finish. The scalable form is a fitting alignment: # one pass that lets the match begin anywhere in t for free, so every end # position is scored at once instead of every pair being scored separately. let matches = range(0, len(sequence)) |> flat_map(|start| range(1, len(sequence) - start + 1) |> filter(|length| edit_distance(motif, substr(sequence, start, length)) <= k) |> map(|length| { start: start + 1, length: length })) println("Result:") for hit in matches { println(" " + str(hit.start) + " " + str(hit.length)) } println("Expected: 1 4 / 1 5 / 1 6") fn test_ksim_finding_all_similar_motifs() { let written = matches |> map(|m| str(m.start) + " " + str(m.length)) assert sort(written) == sort(["1 4", "1 5", "1 6"]), "KSIM: got " + join(written, " / ") # Every reported substring really is within k, and every one omitted is not — # the second half is the part a filter can silently get wrong. for hit in matches { let piece = substr(sequence, hit.start - 1, hit.length) assert edit_distance(motif, piece) <= k, "KSIM: " + piece + " is further than " + str(k) } let missed = range(0, len(sequence)) |> flat_map(|start| range(1, len(sequence) - start + 1) |> filter(|length| edit_distance(motif, substr(sequence, start, length)) <= k) |> filter(|length| (matches |> count_if(|m| m.start == start + 1 and m.length == length)) == 0)) assert len(missed) == 0, "KSIM: a qualifying substring was not reported" # An exact search finds nothing here, which is the whole point of allowing k. assert contains(sequence, motif) == false, "KSIM: the motif does not occur exactly" assert len(matches) > 0, "KSIM: yet three approximate matches exist" } ``` ## RSUB — Identifying Reversing Substitutions [Problem statement](https://rosalind.info/problems/rsub/) ```biolang # Rosalind: RSUB — Identifying Reversing Substitutions # https://rosalind.info/problems/rsub/ # # Given: A rooted binary tree with every node labelled by a string. # Return: Every reversing substitution — a change that later changes back, with # no other change in between. # (((ostrich,cat)rat,mouse)dog,elephant)robot; let children = { "robot": ["dog", "elephant"], "dog": ["rat", "mouse"], "rat": ["ostrich", "cat"], } let sequences = { "robot": "AATTG", "dog": "GGGCA", "mouse": "AAGAC", "rat": "GTTGT", "cat": "GAGGC", "ostrich": "GTGTC", "elephant": "AATTC", } let root = "robot" # A site that mutates and later mutates back looks, from the tips alone, as # though nothing ever happened — the ancestor and the descendant agree. Only the # internal labels reveal it, which is why this problem hands them over rather # than asking for them. # # It matters because reversing substitutions are exactly what makes distant # relationships hard to recover: the signal erases itself, and two lineages look # more similar than their history warrants. fn kids_of(node) { if has_key(children, node) then children[node] else [] } # Walk down from a node where a change occurred, following only lineages that # still carry the new character, and report wherever it changes back. fn reversions_below(start, position, was, became) { let found = [] let frontier = kids_of(start) while len(frontier) > 0 { let node = frontier[0] frontier = slice(frontier, 1, len(frontier)) let here = substr(sequences[node], position, 1) if here == was { # Changed back, with nothing else in between. found = push(found, node) } else { # Still carrying the substitution; keep descending. Any third # character ends the lineage's relevance. if here == became { frontier = frontier + kids_of(node) } } } found } let width = len(sequences[root]) let reported = [] let stack = [root] while len(stack) > 0 { let parent = stack[len(stack) - 1] stack = slice(stack, 0, len(stack) - 1) for child in kids_of(parent) { stack = push(stack, child) for i in range(0, width) { let was = substr(sequences[parent], i, 1) let became = substr(sequences[child], i, 1) if was != became { for reverted in reversions_below(child, i, was, became) { reported = push(reported, child + " " + reverted + " " + str(i + 1) + " " + was + "->" + became + "->" + was) } } } } } println("Result:") for line in sort(reported) { println(" " + line) } println("Expected (any order): dog mouse 1 A->G->A / dog mouse 2 A->G->A") println(" rat ostrich 3 G->T->G / rat cat 3 G->T->G / dog rat 3 T->G->T") fn test_rsub_reversing_substitutions() { let expected = ["dog mouse 1 A->G->A", "dog mouse 2 A->G->A", "rat ostrich 3 G->T->G", "rat cat 3 G->T->G", "dog rat 3 T->G->T"] assert sort(reported) == sort(expected), "RSUB: got " + join(sort(reported), " / ") # Every report must describe a real reversion: the character genuinely # differs from the parent and genuinely returns at the named descendant. for line in reported { let parts = split(line, " ") let changed = parts[0] let reverted = parts[1] let position = int(parts[2]) - 1 assert substr(sequences[changed], position, 1) != substr(sequences[reverted], position, 1), "RSUB: " + line + " does not actually revert" } # The tips alone hide this: robot and mouse agree at position 1 despite two # substitutions on the path between them. assert substr(sequences["robot"], 0, 1) == substr(sequences["mouse"], 0, 1), "RSUB: the endpoints agree, which is what makes the change invisible" assert substr(sequences["dog"], 0, 1) != substr(sequences["robot"], 0, 1), "RSUB: yet the intermediate differs from both" } ``` ## REAR — Reversal Distance [Problem statement](https://rosalind.info/problems/rear/) ```biolang # Rosalind: REAR — Reversal Distance # https://rosalind.info/problems/rear/ # # Given: Up to five pairs of permutations of length 10. # Return: The reversal distance for each pair. let pairs = [ { from_order: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], to_order: [3, 1, 5, 2, 7, 4, 9, 6, 10, 8] }, { from_order: [3, 10, 8, 2, 5, 4, 7, 1, 6, 9], to_order: [5, 2, 3, 1, 7, 4, 10, 8, 6, 9] }, { from_order: [8, 6, 7, 9, 4, 1, 3, 10, 2, 5], to_order: [8, 2, 7, 6, 9, 1, 5, 3, 10, 4] }, { from_order: [3, 9, 10, 4, 1, 8, 6, 7, 5, 2], to_order: [2, 9, 8, 5, 1, 7, 3, 4, 6, 10] }, { from_order: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], to_order: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] }, ] # The 2-break distance in BA6C had a closed form — blocks minus cycles. Reversal # distance has no such formula at this size, so it has to be searched for, and # the search is why `reversal_distance` is a builtin rather than written here: # ten elements admit 45 reversals and 3.6 million reachable orders, and the # distance reaches 9. A one-sided search to that depth is hopeless; the builtin # searches from both ends and meets in the middle, so each side only reaches # depth 4 or 5. let distances = pairs |> map(|pair| reversal_distance(pair.from_order, pair.to_order)) println("Result: " + (distances |> map(|d| str(d)) |> join(" "))) println("Expected: 9 4 5 7 0") fn test_rear_reversal_distance() { assert (distances |> map(|d| str(d)) |> join(" ")) == "9 4 5 7 0", "REAR: got " + str(distances) # Identical permutations are no distance apart, and the measure is symmetric # because a reversal undoes itself. assert distances[4] == 0, "REAR: the last pair is already sorted" for pair in pairs { assert reversal_distance(pair.to_order, pair.from_order) == reversal_distance(pair.from_order, pair.to_order), "REAR: the distance must be symmetric" } } ``` ## SORT — Sorting by Reversals [Problem statement](https://rosalind.info/problems/sort/) ```biolang # Rosalind: SORT — Sorting by Reversals # https://rosalind.info/problems/sort/ # # Given: Two permutations of length 10. # Return: The reversal distance, and a shortest collection of reversals taking # the first to the second. let from_order = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] let to_order = [1, 8, 9, 3, 2, 7, 6, 5, 4, 10] # REAR asks how far apart two gene orders are; this asks for the route. The # search is the same — the distance is just the length of what comes back — but # the sequence is what actually says how the rearrangement happened, and a # distance without one is unfalsifiable. # # Any shortest sequence is correct. Reversals commute in some orders and not # others, so several routes of the same length usually exist. let steps = sorting_reversals(from_order, to_order) println("Result: " + str(len(steps))) for step in steps { println(" " + str(step[0]) + " " + str(step[1])) } println("Expected: 2, then 4 9 and 2 5 — any shortest route is accepted") fn test_sort_sorting_by_reversals() { assert len(steps) == 2, "SORT: expected 2 reversals, got " + str(len(steps)) assert len(steps) == reversal_distance(from_order, to_order), "SORT: the route must be as long as the distance" # The route has to work. Applying the reversals in order must produce the # target — a plausible-looking list that does not sort is the failure this # problem invites. let current = from_order for step in steps { let from_index = step[0] - 1 let to_index = step[1] - 1 let head = range(0, from_index) |> map(|i| current[i]) let middle = range(from_index, to_index + 1) |> map(|i| current[to_index - (i - from_index)]) let tail = range(to_index + 1, len(current)) |> map(|i| current[i]) current = head + middle + tail } assert current == to_order, "SORT: the reversals give " + str(current) + ", not " + str(to_order) # Every interval is a real one, and reverses more than a single element. for step in steps { assert step[0] >= 1 and step[1] <= len(from_order), "SORT: interval out of range" assert step[0] < step[1], "SORT: reversing one element does nothing" } } ```