Bioinformatics Algorithms — Problems &

Solutions Guide
Chapter 2 · Chapter 4 · Chapter 9 (Lecture Slides)
Table of Contents
1. Chapter 2 — Problems & Solutions
2. Chapter 4 — Problems & Solutions
3. Chapter 9 — Lecture Summary & Worked Examples

Chapter 2 — Algorithms & Complexity

{#chapter-2}
Problem 2.1 — Min and Max in a List
Problem: Write an algorithm that, given a list of n numbers, returns the largest and smallest numbers.
Can you do it in only 3n/2 comparisons?
Solution:
The standard approach scans twice → 2(n−1) comparisons. The optimal approach processes pairs of
elements:
 FIND_MIN_MAX(A, n):
if n is odd:
min = max = A[0]
start = 1
else:
if A[0] < A[1]: min = A[0]; max = A[1]
else: min = A[1]; max = A[0]
start = 2

for i = start to n−1 step 2:

if A[i] < A[i+1]:
if A[i] < min: min = A[i]
if A[i+1] > max: max = A[i+1]
else:
if A[i+1] < min: min = A[i+1]
if A[i] > max: max = A[i]

return (min, max)

Analysis:
Each pair uses exactly 3 comparisons (1 internal + 1 vs min + 1 vs max)
n/2 pairs → 3n/2 total comparisons ✓
Time complexity: O(n)

Problem 2.2 — Iterating Over Multi-Dimensional Index

Problem: Write recursive and iterative algorithms to iterate over every index from (0,0,…,0) to
(n₁,n₂,…,nₐ).
Solution — Iterative (Odometer style):
 ITERATE_ITERATIVE(n[1..d]):
index[1..d] = [0, 0, ..., 0]
while True:
process(index)
pos = d
while pos >= 1:
index[pos] = index[pos] + 1
if index[pos] <= n[pos]: break
index[pos] = 0
pos = pos − 1
if pos == 0: break

Solution — Recursive:
ITERATE_RECURSIVE(n[1..d], index[1..d], dim):
if dim > d:
process(index)
return
for i = 0 to n[dim]:
index[dim] = i
ITERATE_RECURSIVE(n, index, dim + 1)

// Call as: ITERATE_RECURSIVE(n, [0..0], 1)

Problem 2.3 — Big-O/Ω/Θ of log n vs n

Problem: Is log n = O(n)? Is log n = Ω(n)? Is log n = Θ(n)?
Solution:
Relationship Answer Reason
log n = O(n) YES log n ≤ n for all n ≥ 1
log n = Ω(n) NO log n / n → 0 as n → ∞
log n = Θ(n) NO Requires both O(n) and Ω(n); Ω fails
Conclusion: log n = O(n) only — log n grows strictly slower than n.
Problem 2.4 — Find the Missing Integer
Problem: Given an unsorted list of n−1 distinct integers from range 1 to n, find the missing integer in
linear time.
Solution:
The sum of 1 to n is n(n+1)/2. Subtract the actual sum to get the missing number:
FIND_MISSING(A[1..n−1], n):
expected_sum = n × (n + 1) / 2
actual_sum = 0
for i = 1 to n−1:
actual_sum = actual_sum + A[i]
return expected_sum − actual_sum

Time complexity: O(n) — a single pass through the array.

Problem 2.5 — Space-Efficient Fibonacci

Problem: Modify FIBONACCI(n) to use O(1) space instead of O(n).
Solution:
Only the two most recent values are ever needed:
FIBONACCI_CONSTANT_SPACE(n):
if n == 0: return 0
if n == 1: return 1
prev = 0
curr = 1
for i = 2 to n:
next = prev + curr
prev = curr
curr = next
return curr

Space: O(1)
Time: O(n)

Problem 2.6 — Closed-Form Fibonacci (Binet's Formula)

Problem: Prove that Fₙ = (1/√5)(φⁿ − φ̄ⁿ) where φ = (1+√5)/2, φ̄ = (1−√5)/2.
Proof by Strong Induction:
Base cases:
F₀ = (1/√5)(1 − 1) = 0 ✓
F₁ = (1/√5)(φ − φ̄) = (1/√5)(√5) = 1 ✓
Inductive step: Assume true for all k ≤ n. Both φ and φ̄ satisfy x² = x + 1, so φ² = φ+1 and φ̄² = φ̄+1.
Fₙ₊₁ = Fₙ + Fₙ₋₁
= (1/√5)(φⁿ − φ̄
ⁿ) + (1/√5)(φⁿ⁻¹ − φ̄
ⁿ⁻¹)
= (1/√5)[φⁿ⁻¹(φ+1) − φ̄
ⁿ⁻¹(φ̄
+1)]
= (1/√5)[φⁿ⁻¹·φ² − φ̄
ⁿ⁻¹·φ̄
²] ← since φ+1 = φ²
= (1/√5)(φⁿ⁺¹ − φ̄
ⁿ⁺¹) ✓

Problem 2.7 — Fibonacci Faster than O(n)

Problem: Design an algorithm to compute the nth Fibonacci number in less than O(n) time.
Solution — Matrix Exponentiation:
Using the identity: [[1,1],[1,0]]ⁿ = [[Fₙ₊₁, Fₙ],[Fₙ, Fₙ₋₁]]
MATRIX_POWER(M, n):
if n == 1: return M
if n is even:
half = MATRIX_POWER(M, n/2)
return MATRIX_MULTIPLY(half, half)
else:
return MATRIX_MULTIPLY(M, MATRIX_POWER(M, n−1))

FAST_FIBONACCI(n):
if n == 0: return 0
M = [[1,1],[1,0]]
result = MATRIX_POWER(M, n)
return result[0][1]

Time complexity: O(log n) matrix multiplications, each O(1) for 2×2 matrices.
Problem 2.8 — Limited Lifespan Rabbit Model
Problem: Model rabbit population where rabbits die after k years. Write a recurrence and pseudocode.
Recurrence:
R(n) = R(n−1) + R(n−2) − R(n − ⌊k⌋ − 1)

(Add new births, subtract pairs dying at age ⌊k⌋+1)

LIMITED_FIBONACCI(n, k):
R[0] = 1; R[1] = 1
lifespan = floor(k) + 1
for i = 2 to n:
if i <= lifespan:
R[i] = R[i−1] + R[i−2]
else:
R[i] = R[i−1] + R[i−2] − R[i − lifespan − 1]
return R[n]

For k=2.5 and the sequence gives: 1, 1, 2, 3, 4 — growing more slowly than standard Fibonacci.

Problem 2.9 — Iterative Tower of Hanoi

Problem: Write a non-recursive algorithm for Tower of Hanoi.
Solution:
Total moves = 2ⁿ − 1. For odd n, cycle pegs [source, target, auxiliary]; for even n, [source, auxiliary, target].
 ITERATIVE_HANOI(n, source, target, auxiliary):
total_moves = 2^n − 1
if n is odd:
pegOrder = [source, target, auxiliary]
else:
pegOrder = [source, auxiliary, target]

for move = 1 to total_moves:

if move mod 3 == 1:
MOVE_DISK(pegOrder[0], pegOrder[1])
elif move mod 3 == 2:
MOVE_DISK_NON_SMALLEST(pegOrder[0], pegOrder[2])
else:
MOVE_DISK_NON_SMALLEST(pegOrder[1], pegOrder[2])

Time: O(2ⁿ) — optimal, since exactly 2ⁿ−1 moves are required.

Problem 2.10 — Sum of First n Integers

Problem: Prove that Σᵢ₌₁ⁿ i = n(n+1)/2.
Proof by Induction:
Base case: n=1: LHS=1, RHS=1(2)/2=1 ✓
Inductive step: Assume Σᵢ₌₁ᵏ i = k(k+1)/2. Then:
Σᵢ₌₁ᵏ⁺¹ i = k(k+1)/2 + (k+1)
= (k+1)(k/2 + 1)
= (k+1)(k+2)/2 ✓

This is exactly the formula with n = k+1. □

Problem 2.11 — Geometric Series Proofs

Problem: Prove Σᵢ₌₁ⁿ 2ⁱ = 2ⁿ⁺¹ − 2 and Σᵢ₌₁ⁿ 2⁻ⁱ = 1 − 2⁻ⁿ.
Proof:
From the geometric series formula: Σᵢ₌₀ⁿ rⁱ = (rⁿ⁺¹−1)/(r−1)
Part 1 (r=2):
 Σᵢ₌₁ⁿ 2ⁱ = 2·(2ⁿ−1)/(2−1) = 2ⁿ⁺¹ − 2 ✓

Part 2 (r=1/2):
Σᵢ₌₁ⁿ (1/2)ⁱ = (1/2)·(1−(1/2)ⁿ)/(1−1/2) = 1 − (1/2)ⁿ = 1 − 2⁻ⁿ ✓

Problem 2.12 — Fix BetterChange With a New

Denomination
Problem: BETTERCHANGE fails on {25,20,10,5,1}. Add a denomination to make it correct for all M.
Solution:
The failure occurs because 25 is not a multiple of 20. For M=40, greedy picks 25+10+5=3 coins instead of
optimal 20+20=2.
Fix: Add denomination 40 → set becomes {40,25,20,10,5,1}.
Now M=40 → 1 coin directly. The greedy algorithm never needs to "commit" to 25 when 40 is available.
Key insight: A greedy algorithm is correct for a denomination set when each denomination divides
the next larger one, or when the "greedy choice property" can be proven.

Problem 2.13 — Average Coins by UCHANGE(M)

Problem: Design an algorithm that computes the average number of coins returned by UCHANGE(M) as M
varies from 1 to 100.
AVERAGE_COINS(denominations):
total_coins = 0
for M = 1 to 100:
coins = UCHANGE(M, denominations)
total_coins = total_coins + coins
return total_coins / 100

Time: O(100 × d) = O(d) where d = number of denominations.

Problem 2.14 — Verify Correctness of BetterChange

Problem: Given arbitrary denominations, decide if BETTERCHANGE is correct.
VERIFY_BETTERCHANGE(c[1..d]):
SORT_DESC(c)
for M = 1 to c[1] × d:
greedy_count = BETTERCHANGE(M, c)
optimal_count = DP_CHANGE(M, c)
if greedy_count != optimal_count:
return FALSE
return TRUE

Test up to the LCM of denominations — if greedy fails for any M in that range, it's incorrect.

Problem 2.15 — King on a Chessboard Game

Problem: King moves right, down, or diagonally SE on 8×8 board. Two players alternate; player who
places king on lower-right wins. Who wins?
Solution:
Label positions WINNING (W) or LOSING (L) by working backwards:
L = ALL reachable positions are W (opponent always wins from there)
W = at least ONE reachable position is L
Pattern: Position (i,j) is LOSING iff i ≡ j (mod 2) — equal parity.
Position (1,1): both odd → same parity → LOSING for Player 1.
Conclusion: The SECOND player wins. After any Player 1 move, parity changes; Player 2 restores
equal parity, eventually landing on (8,8).

Problem 2.16 — Rock Splitting Game

Problem: n rocks, one pile. Players alternate splitting any pile of >1 rocks. Last player to move wins. Who
wins?
Solution:
Each split increases the total number of piles by exactly 1. Starting from 1 pile, exactly n−1 total splits
occur before all piles have size 1.
n even → n−1 is odd → Player 1 makes the last move → Player 1 wins
 n odd → n−1 is even → Player 2 makes the last move → Player 2 wins

Problem 2.17 — Virus vs. Bacteria

Problem: n bacteria, 1 virus. Each minute: viruses kill 1 bacterium each and double; surviving bacteria
also double. Will viruses win?
Solution:
Let B(t) = bacteria, V(t) = viruses at time t.
V(t) = 2^t
B(t) = 2^t × (n − t)

B(t) = 0 when t = n. Yes, viruses kill all bacteria after exactly n minutes.
VIRUS_KILLS(n):
V = 1; B = n; t = 0
while B > 0:
V_next = 2 × V
B_next = 2 × (B − V)
V = V_next; B = B_next; t = t + 1
return t // returns n

Time: O(n) — exactly n iterations.

Problem 2.18 — Identifying Honest Professors

Problem: 100 professors (>50 honest). Find all honest ones with ≤198 questions ("Is X honest?" asked to
Y).
Solution — Tournament Algorithm:
Key insight: If A says B is honest AND B says A is honest → they are the same type (both honest or both
dishonest).
Phase 1 (100 questions): Pair up professors (50 pairs). Ask each pair about each other (2 × 50 = 100
questions). Discard any pair where either gives a "no" answer. Keep pairs where both say "yes" — take one
representative from each surviving pair.
Phase 2 (≤98 questions): The majority of survivors are still honest. Ask all about one candidate; majority
vote reveals truth. Classify all 100 professors.
Total: ≤ 198 questions. □
Problem 2.19 — Transform 8×8 Table to All Zeros
Problem: Operations: double all values in a row, or subtract 1 from all values in a column. Reach all zeros.
Solution:
Strategy: For each column j:
1. Find the maximum value M in that column
2. For each row i, double row i until T[i][j] reaches M (power-of-2 alignment)
3. Subtract column j exactly M times
Since all entries are natural numbers and we can always double to match a power of 2, the algorithm
terminates.
Running time: O(max_value × 8²) column operations.

Problem 2.20 — Chameleon Color Convergence

Problem: n black, m green, k brown chameleons. When two different-colored meet, both become the
third color. Can all become one color?
Solution:
Invariant: (n−m) mod 3, (m−k) mod 3, and (n−k) mod 3 are all preserved under each meeting event.
Condition for convergence: All same color is possible if and only if n ≡ m ≡ k (mod 3).
Example n=1, m=3, k=5: 1 mod 3=1, 3 mod 3=0, 5 mod 3=2. Not all equal → impossible.

Chapter 4 — Exhaustive Search {#chapter-4}

Problem 4.1 — Compute Delta(X)
Problem: Write an algorithm that, given a set X, calculates the multiset ΔX.
Definition: ΔX = {|xᵢ − xⱼ| : xᵢ, xⱼ ∈ X, i ≠ j}
 COMPUTE_DELTA(X[1..n]):
delta = empty multiset
for i = 1 to n:
for j = 1 to n:
if i != j:
[Link](|X[i] − X[j]|)
return delta

// Time: O(n²)
// If X has n elements → ΔX has n(n−1) elements (ordered)

Problem 4.2 — Solve the Partial Digest Problem

Problem: Given L = {1,1,1,2,2,3,3,3,4,4,5,5,6,6,6,9,9,10,11,12,15}, find X such that ΔX = L.
Solution (using the PARTIALDGEST algorithm):
width = max(L) = 15. So X must contain {0, …, 15}.
Step max(L) Try y Δ(y, X) Valid? Action
1 15 width=15 — — X={0,15}, remove 15
2 12 y=12 {12,3} ✓ both in L X={0,12,15}, remove {12,3}
3 11 y=11 {11,1,4} ✓ X={0,11,12,15}, remove {1,4,11}
4 10 y=10 {10,1,2,5} ✓ X={0,10,11,12,15}, remove {1,2,5,10}
5 9 y=9 {9,1,2,3,6} ✓ X={0,9,10,11,12,15}, remove {1,2,3,6,9}
6 6 y=6 {6,3,4,5,6,9} ✓ X={0,6,9,10,11,12,15} — DONE
Answer: X = {0, 6, 9, 10, 11, 12, 15} (or its reflection {0, 3, 4, 5, 6, 9, 15})

Problem 4.3 — Generate All m-element Subsets

Problem: Write an algorithm that generates all m-element subsets of an n-element set.
 SUBSETS(S[1..n], m):
GENERATE(S, m, 1, [])

GENERATE(S, m, start, current):

if length(current) == m:
OUTPUT(current)
return
if start > n: return

// Include S[start]
[Link](S[start])
GENERATE(S, m, start+1, current)
[Link]()

// Exclude S[start]
GENERATE(S, m, start+1, current)

Time: O(C(n,m) × m) where C(n,m) = n!/(m!(n−m)!) — optimal since every subset must be visited.

Problem 4.4 — m-element Subsets of a Multiset

Problem: Generate all m-element subsets of an n-element multiset (with repeated elements).
MULTISET_SUBSETS(S[1..n], m):
SORT(S)
GENERATE_MULTI(S, m, 1, [])

GENERATE_MULTI(S, m, start, current):

if length(current) == m:
OUTPUT(current)
return
prev = NULL
for i = start to n − (m − length(current)) + 1:
if S[i] != prev: // skip duplicates at same recursion level
[Link](S[i])
GENERATE_MULTI(S, m, i+1, current)
[Link]()
prev = S[i]

Example {1,2,2,3}, m=2: Produces {1,2}, {1,3}, {2,2}, {2,3} — 4 subsets, no duplicates.
Problem 4.5 — Homometric Sets (Proof)
Problem: Prove that U⊕V = {u+v} and U⊖V = {u−v} are homometric for any U, V.
Proof:
Two sets A and B are homometric if ΔA = ΔB (same multiset of pairwise differences).
Let a₁=u₁+v₁, a₂=u₂+v₂ ∈ U⊕V. Then a₁−a₂ = (u₁−u₂)+(v₁−v₂).
Let b₁=u₁−v₁, b₂=u₂−v₂ ∈ U⊖V. Then b₁−b₂ = (u₁−u₂)−(v₁−v₂).
The multiset {(u₁−u₂)+(v₁−v₂)} over all pairs equals the multiset {(u₁−u₂)−(v₁−v₂)} because for every (u,v)
difference pair, both +v and −v contributions appear symmetrically in ΔV.
Therefore Δ(U⊕V) = Δ(U⊖V), proving homometricity. □

Problem 4.6 — Generating Functions for Delta

Problem: Verify that ΔA(x) = A(x)·A(x⁻¹) where A(x) = Σ x^(aᵢ).
Solution:
A(x)·A(x⁻¹) = (Σᵢ x^aᵢ)(Σⱼ x^(−aⱼ)) = Σᵢ Σⱼ x^(aᵢ−aⱼ)

This is exactly the generating function of ΔA = {aᵢ−aⱼ : i≠j} — confirmed. ✓

For U⊕V and U⊖V:
Δ(U⊕V)(x) = U(x)V(x) · U(x⁻¹)V(x⁻¹) = [U(x)U(x⁻¹)][V(x)V(x⁻¹)]
Δ(U⊖V)(x) = U(x)V(x⁻¹) · U(x⁻¹)V(x) = [U(x)U(x⁻¹)][V(x⁻¹)V(x)]

They are equal → confirming homometricity. ✓

Problem 4.7 — Compact PARTIALDGEST Pseudocode

Problem: Write pseudocode for PARTIALDGEST in fewer lines.
 PARTIAL_DIGEST(L):
width = max(L)
[Link](width)
X = {0, width}
PLACE(L, X)

PLACE(L, X):
if L is empty: OUTPUT(X); return
y = max(L)
if VALID(y, X, L):
[Link](y); L.remove_all(delta(y, X))
PLACE(L, X)
[Link](y); L.restore_all(delta(y, X))
y2 = width − y
if VALID(y2, X, L):
[Link](y2); L.remove_all(delta(y2, X))
PLACE(L, X)
[Link](y2); L.restore_all(delta(y2, X))

VALID(y, X, L):
for each x in X:
if |y − x| not in L: return FALSE
return TRUE

Problem 4.8 — Smallest Non-Unique Delta(X)

Problem: Find a ΔX with the smallest number of elements that could arise from more than one X.
Solution:
The smallest such example requires 4 points. All sets with 2 or 3 points produce unique ΔX (up to
reflection and shift).
Example of a homometric pair:
X₁ = {0, 1, 5, 11, 13}
X₂ = {0, 1, 8, 9, 13}
Both produce the same multiset ΔX despite being different sets. The minimum non-trivial case has 4
points with ΔX having 6 elements.

Problem 4.9 — Double Digest Problem (DDP)

Problem: Devise brute force and branch-and-bound algorithms for DDP.
Brute Force:
BRUTE_FORCE_DDP(A, B, AB):
for each permutation P_A of A:
for each permutation P_B of B:
if CONSISTENT(P_A, P_B, AB):
return (P_A, P_B)

Complexity: O((|A|!)² ) — very slow.

Branch-and-Bound Improvement:
BRANCH_BOUND_DDP(A, B, AB):
build restriction map incrementally:
for each cut site of enzyme A (left to right):
place next site at position p
compute partial double-digest fragments
if any required fragment is absent from AB:
PRUNE this branch and backtrack
else:
continue placing sites

Key: Prune whenever the partial assignment creates a fragment length not present in AB.

Problem 4.10 — Probed Partial Digest Problem (PPDP)

Problem: Design brute force and branch-and-bound algorithms for PPDP.
Setup: Probe at position 0. A = negative sites, B = positive sites. Input: {b−a : a∈A, b∈B}.
Branch-and-Bound:
 EXTEND(A, B, Remaining):
if Remaining is empty: OUTPUT(A, B); return
for candidate in CANDIDATE_SITES(Remaining):
if VALID_POSITIVE(candidate, A, B, Remaining):
[Link](candidate)
new_remaining = Remaining − {b−a : a in A}
EXTEND(A, B, new_remaining)
[Link](candidate)
if VALID_NEGATIVE(candidate, A, B, Remaining):
[Link](candidate)
new_remaining = Remaining − {b−a : b in B}
EXTEND(A, B, new_remaining)
[Link](candidate)

Problem 4.11 — Complete k-ary Tree Vertex Count

Problem: Find a closed-form expression for total vertices in a complete balanced k-ary tree of height L.
Solution:
Depth Nodes
0 (root) 1
1 k
2 k²
… …
L (leaves) k^L
L k L+1 −1
Total = ∑d=0 k d =​

k−1
​ (k 
= 1)

Example: DNA search tree (k=4) of height 8:

Total = (4⁹ − 1)/3 = 262143/3 = 87,381 nodes

Problem 4.12 — String Pattern Search

Problem: Given text T and pattern s, find the first occurrence of s in T.
 FIND_PATTERN(T, s):
n = length(T); m = length(s)
for i = 1 to n − m + 1:
match = TRUE
for j = 1 to m:
if T[i+j−1] != s[j]:
match = FALSE; break
if match: return i
return −1

Naive complexity: O(n × m)

Improved: KMP algorithm achieves O(n + m)

Problem 4.13 — Approximate Pattern Search (Hamming

Distance)
Problem: Find first occurrence of s' in T with Hamming distance dH(s, s') ≤ k.
APPROX_FIND(T, s, k):
n = length(T); m = length(s)
for i = 1 to n − m + 1:
hamming = 0
for j = 1 to m:
if T[i+j−1] != s[j]:
hamming = hamming + 1
if hamming > k: break // early exit pruning
if hamming <= k: return i
return −1

Complexity: O(n × m) worst case. The early break significantly improves average case.

Problem 4.14 — Count l-mer Frequencies

Problem: Count occurrences of each l-mer in a string of length n.
 COUNT_LMERS(text, l):
freq = empty hash map
n = length(text)
for i = 1 to n − l + 1:
lmer = text[i .. i+l−1]
freq[lmer] = [Link](lmer, 0) + 1
return freq

Complexity: O(n × l) using hash map (O(l) per hash).

For a bacterial genome (~4×10⁶ bp) with l=6: ~4M iterations. Random strings show uniform
distribution (~n/4^l per l-mer); real genomes show significant deviation — over/under-represented
motifs are biologically significant.

Problem 4.15 — Identify ANOTHERMOTIFSEARCH

Problem: Identify which algorithm ANOTHERMOTIFSEARCH is a cousin of.
Answer: ANOTHERMOTIFSEARCH is a cousin of Exhaustive Brute-Force Motif Search.
ANOTHERMOTIFSEARCH(DNA, t, n, l):
s ← (1, 1, ..., 1)
bestMotif ← FINDINSEQ(s, 1, t, n, l)
return bestMotif

FINDINSEQ(s, currentSeq, t, n, l):

bestScore ← 0
for j ← 1 to n − l + 1:
s[currentSeq] ← j
if currentSeq != t:
s ← FINDINSEQ(s, currentSeq + 1, t, n, l)
if Score(s) > bestScore:
bestScore ← Score(s)
bestMotif ← s
return bestMotif

Aspect ANOTHERMOTIFSEARCH Standard Exhaustive

Approach Recursive DFS Iterative nested loops
Explores All (n−l+1)^t combos All (n−l+1)^t combos
 Pruning None None
Result Optimal Optimal
Complexity O((n−l+1)^t × lt) O((n−l+1)^t × lt)

Problem 4.16 — Identify YETANOTHERMOTIFSEARCH

Problem: Identify which algorithm YETANOTHERMOTIFSEARCH is a cousin of.
Answer: YETANOTHERMOTIFSEARCH is a cousin of Branch-and-Bound Motif Search.
FIND(s, currentSeq, t, n, l):
i ← currentSeq
bestScore ← 0
for j ← 1 to n − l + 1:
si ← j
bestPossibleScore ← Score(s, i) + (t − i) × l // UPPER BOUND
if bestPossibleScore > bestScore: // PRUNING CHECK
if currentSeq != t:
s ← FIND(s, currentSeq + 1, t, n, l)
if Score(s) > bestScore:
bestScore ← Score(s)
bestMotif ← s
return bestMotif

Key difference from ANOTHERMOTIFSEARCH:

The bound bestPossibleScore = Score(s,i) + (t−i)×l is an upper bound on the best achievable
score from the current partial assignment. If this bound ≤ current bestScore, the entire subtree is pruned
without exploration.

Problem 4.17 — Tighter Bound for Median String B&B

Problem: Derive a tighter bound using the split of l-mer w into prefix u and suffix v.
Proof:
Split motif w of length l into prefix u (length p) and suffix v (length l−p).
For any sequence j and position i:
dH(w, s_{j,i}) = dH(u, s_{j,i}[1..p]) + dH(v, s_{j,i}[p+1..l])
Therefore:
TotalDistance(w, DNA) = Σⱼ minᵢ dH(w, s_{j,i})
≥ Σⱼ minᵢ [dH(u, s_{j,i}[1..p]) + dH(v, s_{j,i}[p+1..l])]
≥ Σⱼ [minᵢ dH(u, s_{j,i}[1..p]) + minᵢ dH(v, s_{j,i}
[p+1..l])]
= TotalDistance(u, DNA) + TotalDistance(v, DNA)

Tighter bound: If TotalDistance(u, DNA) + TotalDistance(v, DNA) > bestScore , prune the
entire subtree of l-mers beginning with u. Since TotalDistance(u) can be precomputed for all (l/2)-mers,
this prune is cheap to compute and substantially tighter than the basic bound. □

Chapter 9 — Randomized Algorithms & Motif

Finding {#chapter-9}
Summary: Key Concepts
Why Randomized Algorithms?
Motif finding is NP-hard. No fast exact algorithm is known. Randomized algorithms:
1. Make random decisions so no input consistently triggers worst-case behavior
2. Can be run many times — at least one run will almost certainly find the correct answer

The Motif Finding Pipeline

Starting positions s = (s₁,...,sₜ)
↓
Alignment matrix (t × l)
↓
COUNT(Motifs) — 4 × l count matrix
↓ (divide by t)
PROFILE(Motifs) — 4 × l frequency matrix
↓
Consensus string (most frequent per column)
↓
Score(Motifs) = Σⱼ (t − max_b COUNT[b,j])
Profile Scoring
Formula: Prob(a|P) = Π_{j=1}^{l} P[a_j, j]
Worked Example:
Profile P: Pos1 Pos2 Pos3 Pos4 Pos5 Pos6
A: 1/2 7/8 3/8 0 1/8 0
C: 1/8 0 1/2 5/8 3/8 0
T: 1/8 1/8 0 0 1/4 7/8
G: 1/4 0 1/8 3/8 1/4 1/8

Prob(aaacct|P) = 1/2 × 7/8 × 3/8 × 5/8 × 3/8 × 7/8 = 0.0336 ← consensus

Prob(atacag|P) = 1/2 × 1/8 × 3/8 × 5/8 × 1/8 × 1/8 = 0.0016 ← lower

Laplace's Rule of Succession (Pseudocounts)

Problem: Any zero in PROFILE makes Prob(a|P) = 0, destroying all other position information.
Fix: Add 1 to every COUNT entry before computing PROFILE.
Before pseudocounts: After adding 1:
Motifs: T A A C A: 2+1 1+1 1+1 1+1 → 3/8 2/8 2/8 2/8
G T C T C: 0+1 1+1 1+1 1+1 → 1/8 2/8 2/8 2/8
A C T A G: 1+1 1+1 1+1 0+1 → 2/8 2/8 2/8 1/8
A G G T T: 1+1 1+1 1+1 2+1 → 2/8 2/8 2/8 3/8

Denominator = t + 4 = 4 + 4 = 8

Greedy Profile Motif Search

 GREEDYMOTIFSEARCH(Dna, k, t):
BestMotifs ← first k-mer from each string in Dna
for each k-mer Motif in Dna[1]:
Motif_1 ← Motif
for i = 2 to t:
Profile ← PROFILE(Motif_1, ..., Motif_{i-1}) // with pseudocounts
Motif_i ← PROFILE_MOST_PROBABLE(Dna[i], k, Profile)
Motifs ← (Motif_1, ..., Motif_t)
if SCORE(Motifs) < SCORE(BestMotifs):
BestMotifs ← Motifs
return BestMotifs

Complexity: O((n−k+1) × t × k × (n−k+1))

Limitation: Greedy early choices may be suboptimal.

Randomized Motif Search

RANDOMIZEDMOTIFSEARCH(Dna, k, t):
Motifs ← randomly select one k-mer from each string
BestMotifs ← Motifs
while True:
Profile ← PROFILE(Motifs) // with pseudocounts
Motifs ← MOTIFS(Profile, Dna) // best k-mer per string
if SCORE(Motifs) < SCORE(BestMotifs):
BestMotifs ← Motifs
else:
return BestMotifs

Run ~1000 times and return the best result across all runs.
Why it works: Random starts create statistical bias toward the true motif. Each run has a significant
probability of capturing at least one implanted motif, which pulls the profile toward the correct solution.

Gibbs Sampler
 GIBBSSAMPLER(Dna, k, t, N):
Motifs ← randomly select one k-mer from each string
BestMotifs ← Motifs
for j = 1 to N:
i ← RANDOM(1, t) // remove one random sequence
Profile ← PROFILE(all Motifs except Motif_i) // from t−1 strings
Motif_i ← PROFILE_RANDOMLY_GENERATED(Profile, Dna[i], k)
// ^ sample from distribution, NOT take maximum
if SCORE(Motifs) < SCORE(BestMotifs):
BestMotifs ← Motifs
return BestMotifs

Critical difference: PROFILE_RANDOMLY_GENERATED samples from the probability distribution — it does

not always pick the maximum probability k-mer. This allows escaping local optima.

Gibbs Sampler — Full Step-by-Step Trace

Input: t=5 sequences, l=8. Starting positions: s=(7,11,9,4,1)
Seq 1 (s=7): GTAAACAA[TATTTAT]AGC
Seq 2 (s=11): AAAATTTACC[TTAGAAG]G
Seq 3 (s=9): CCGTACTGT[CAAGCGTG]G
Seq 4 (s=4): TGA[GTAAACGA]CGTCCCA
Seq 5 (s=1): [TACTTAAC]ACCCTGTCAA

Step 1: Randomly remove Sequence 2.

Step 2: Build Profile P from seqs 1, 3, 4, 5:
Pos: 1 2 3 4 5 6 7 8
A: 1/4 2/4 2/4 3/4 1/4 1/4 1/4 2/4
C: 0 1/4 1/4 0 0 2/4 0 1/4
T: 2/4 1/4 1/4 1/4 2/4 1/4 1/4 1/4
G: 1/4 0 0 0 1/4 0 3/4 0
Consensus: T A A A T C G A

Step 3: Compute Prob(each 8-mer in Seq 2 | P):

 Pos 1 (AAAATTTA): 0.000732 ← highest
Pos 2 (AAATTTAC): 0.000122
Pos 8 (ACCTTAGA): 0.000183
All others: 0

Step 4: Normalize to probability distribution:

Ratios: 6 : 1 : 1.5 (divide by minimum 0.000122)
Sum: 8.5

P(pos 1) = 6/8.5 = 0.706

P(pos 2) = 1/8.5 = 0.118
P(pos 8) = 1.5/8.5 = 0.176

Step 5: Sample from distribution → most likely pos 1 → update s₂=1.

Step 6: Repeat until no improvement.

Algorithm Comparison Table

Property Greedy Randomized Gibbs Sampler
Updates per iteration All t motifs (greedy) All t motifs ONE motif (random)
Selection method Deterministic max Deterministic max Probabilistic sample
Local optima risk High Medium Lower
Speed per run Fastest Fast Slower
Reruns needed No Yes (~1000) Yes (~20 long)
Solution quality Suboptimal Good Best
Escapes bad starts No Somewhat Better

Key Formulas Quick Reference

Formula Description
Prob(a|P) = Π_{j=1}^{l} P[a_j, j] Profile probability of k-mer a
 Score(Motifs) = Σⱼ (t − max_b COUNT[b,j]) Total mismatches from consensus
TotalDist(v,DNA) = Σⱼ minᵢ dH(v, s_{j,i}) Median string objective
Nodes = (k^{L+1} − 1) / (k−1) Vertices in complete k-ary tree of height L
Hanoi(n) = 2^n − 1 Moves for n-disk Tower of Hanoi
Σᵢ₌₁ⁿ i = n(n+1)/2 Sum of first n integers
Σᵢ₌₀ⁿ rⁱ = (r^{n+1}−1)/(r−1) Geometric series
dH(u,v) = |{j : u_j ≠ v_j}| Hamming distance
|ΔX| = n(n−1)/2 Number of pairwise distances for n-point set

End of Problems & Solutions Guide