Graphs Dynamic Programming

Graph Algorithms & Dynamic Programming

When to Use

• Implementing Needleman-Wunsch (global) or Smith-Waterman (local) pairwise sequence alignment from scratch instead of calling Bio.Align.PairwiseAligner
• Computing shortest/most-reliable paths on a protein-protein interaction or signaling network (BFS k-hop, Dijkstra with log-transformed edge confidences)
• Building a quick single-linkage clustering or simplified phylogenetic tree from a pairwise distance matrix (Kruskal/Prim MST)
• Ordering a gene regulatory cascade, dependency graph, or pipeline of tasks with precedence constraints (topological sort)
• Selecting an optimal subset under a budget/capacity constraint -- gene panel design, primer set selection (0/1 or unbounded knapsack)

Version Compatibility

Pure Python stdlib (collections, heapq) + NumPy >=1.24. No version-sensitive APIs; works on Python >=3.9. For production alignment/graph work, prefer vetted libraries (Bio.Align, networkx, scipy.sparse.csgraph) -- use this skill's code when you need a from-scratch, dependency-light implementation or must understand/modify the recurrence.

Prerequisites

• pip install numpy (only needed for the DP matrices in NW/SW; graph algorithms are pure stdlib)
• Familiarity with Big-O notation and recursion
• For applying results to real sequences/networks, see biopython (production alignment) and networkx (production graph algorithms)

Complexity Reference

Graph representations: Adjacency matrix O(V^2) space for dense + frequent lookups. Adjacency list O(V+E) for sparse (bio networks: PPI ~20k nodes, ~300k edges -- list saves ~170x memory). Edge list for MST/I/O.

DP Framework

1. State: What does dp[i][j] represent?
2. Base case: Smallest subproblem with known answer
3. Recurrence: How does dp[i] depend on smaller states?
4. Answer extraction: Corner, max cell, or traceback?

Memoization vs tabulation: Memoization (top-down, @lru_cache) computes only needed subproblems but risks stack overflow. Tabulation (bottom-up) enables rolling-array space optimization.

Sequence Alignment Notes

• NW (global): Init borders with gap penalties; traceback from bottom-right
• SW (local): Floor cells at 0; traceback from max cell (stops at 0)
• Affine gaps: Three matrices M/X/Y; open cost once, extend per extra base
• Log-space for max-product paths: min(-log(conf)) via Dijkstra

Code Templates

Adjacency List + BFS + DFS

from collections import defaultdict, deque

graph = defaultdict(list)  # {node: [(neighbor, weight), ...]}
graph[u].append((v, w)); graph[v].append((u, w))

def bfs(graph, start):
    dist = {start: 0}
    queue = deque([start])
    while queue:
        u = queue.popleft()
        for v in graph[u]:
            if v not in dist:
                dist[v] = dist[u] + 1
                queue.append(v)
    return dist

def dfs(graph, start):
    visited, stack, order = set(), [start], []
    while stack:
        u = stack.pop()
        if u not in visited:
            visited.add(u); order.append(u)
            stack.extend(v for v in graph[u] if v not in visited)
    return order

Dijkstra

import heapq
def dijkstra(graph, start):
    dist, heap, visited, pred = {start: 0}, [(0, start)], set(), {start: None}
    while heap:
        d, u = heapq.heappop(heap)
        if u in visited: continue
        visited.add(u)
        for v, w in graph[u]:
            if d + w < dist.get(v, float('inf')):
                dist[v] = d + w; pred[v] = u
                heapq.heappush(heap, (dist[v], v))
    return dist, pred

Union-Find + Kruskal MST

class UF:
    def __init__(self, n): self.p = list(range(n)); self.rank = [0]*n
    def find(self, x):
        if self.p[x] != x: self.p[x] = self.find(self.p[x])
        return self.p[x]
    def union(self, x, y):
        px, py = self.find(x), self.find(y)
        if px == py: return False
        if self.rank[px] < self.rank[py]: px, py = py, px
        self.p[py] = px
        if self.rank[px] == self.rank[py]: self.rank[px] += 1
        return True

def kruskal(vertices, edges):
    idx = {v: i for i, v in enumerate(vertices)}
    uf = UF(len(vertices)); mst, total = [], 0
    for u, v, w in sorted(edges, key=lambda e: e[2]):
        if uf.union(idx[u], idx[v]):
            mst.append((u, v, w)); total += w
            if len(mst) == len(vertices) - 1: break
    return mst, total

Topological Sort (Kahn)

def topo_kahn(vertices, adj):
    in_deg = {v: 0 for v in vertices}
    for u in adj:
        for v in adj[u]: in_deg[v] += 1
    queue = deque(v for v in vertices if in_deg[v] == 0)
    result = []
    while queue:
        u = queue.popleft(); result.append(u)
        for v in adj[u]:
            in_deg[v] -= 1
            if in_deg[v] == 0: queue.append(v)
    return result if len(result) == len(vertices) else None  # None = cycle

0/1 Knapsack (space-optimized + traceback variant)

def knapsack(weights, values, capacity):
    dp = [0] * (capacity + 1)
    for w, v in zip(weights, values):
        for c in range(capacity, w - 1, -1):  # backwards prevents reuse
            dp[c] = max(dp[c], dp[c - w] + v)
    return dp[capacity]

Needleman-Wunsch (global alignment)

import numpy as np
def needleman_wunsch(s1, s2, match=1, mismatch=-1, gap=-2):
    m, n = len(s1), len(s2)
    dp = np.zeros((m+1, n+1), dtype=int)
    dp[:,0] = [i*gap for i in range(m+1)]
    dp[0,:] = [j*gap for j in range(n+1)]
    for i in range(1, m+1):
        for j in range(1, n+1):
            s = match if s1[i-1]==s2[j-1] else mismatch
            dp[i][j] = max(dp[i-1][j-1]+s, dp[i-1][j]+gap, dp[i][j-1]+gap)
    # Traceback
    a1, a2 = [], []; i, j = m, n
    while i>0 or j>0:
        if i>0 and j>0:
            s = match if s1[i-1]==s2[j-1] else mismatch
            if dp[i][j] == dp[i-1][j-1]+s:
                a1.append(s1[i-1]); a2.append(s2[j-1]); i-=1; j-=1; continue
        if i>0 and dp[i][j] == dp[i-1][j]+gap:
            a1.append(s1[i-1]); a2.append('-'); i-=1
        else:
            a1.append('-'); a2.append(s2[j-1]); j-=1
    return dp[m][n], ''.join(reversed(a1)), ''.join(reversed(a2))

Smith-Waterman (local alignment)

def smith_waterman(s1, s2, match=2, mismatch=-1, gap=-1):
    m, n = len(s1), len(s2)
    dp = np.zeros((m+1, n+1), dtype=int)
    best, pos = 0, (0,0)
    for i in range(1, m+1):
        for j in range(1, n+1):
            s = match if s1[i-1]==s2[j-1] else mismatch
            dp[i][j] = max(0, dp[i-1][j-1]+s, dp[i-1][j]+gap, dp[i][j-1]+gap)
            if dp[i][j] > best: best=dp[i][j]; pos=(i,j)
    a1, a2 = [], []; i, j = pos
    while dp[i][j] > 0:
        s = match if s1[i-1]==s2[j-1] else mismatch
        if dp[i][j]==dp[i-1][j-1]+s: a1.append(s1[i-1]); a2.append(s2[j-1]); i-=1; j-=1
        elif dp[i][j]==dp[i-1][j]+gap: a1.append(s1[i-1]); a2.append('-'); i-=1
        else: a1.append('-'); a2.append(s2[j-1]); j-=1
    return best, ''.join(reversed(a1)), ''.join(reversed(a2))

Pitfalls

• Dijkstra with negative weights fails silently -- use Bellman-Ford
• Topo sort on cyclic graph: Kahn returns partial list; check len(result) == len(vertices)
• DFS cycle detection (undirected): track parent to avoid false positives
• Knapsack iterate order: 0/1 backwards; unbounded forwards
• Knapsack traceback requires full 2D table, not space-optimized 1D
• SW traceback stops at 0, not at corner -- do not reuse NW traceback logic
• Affine gap alignment needs 3 matrices (M/X/Y); single matrix gives wrong gap-run scores
• MST on disconnected graph produces a spanning forest, not a single tree

Bioinformatics Applications

See Also

• biopython -- production-grade Bio.Align.PairwiseAligner and BLAST wrappers instead of hand-rolled NW/SW
• networkx -- battle-tested graph library for larger PPI/regulatory networks, shortest paths, and MST
• bio-phylogenetics-distance-calculations -- pairwise distance matrices to feed into MST-based clustering
• bio-gene-regulatory-networks-coexpression-networks -- building the DAGs that topological sort orders