Align the beta- and delta-globin genes
Build affine-gap global and local alignment, follow the traceback, and explain exactly what identity and coverage measure.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the biological question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
RECIPE_ID='02-affine-sequence-alignment'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
Load the embedded snapshot
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'sequences': 'a8ae02a43c841e2932efafaf82a0bf78199baa8921a821fcbd774c8d7649e8b7'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['sequences']
Read the method functions
def affine_align(a, b, local=False, match=2, mismatch=-1, gap_open=-3, gap_extend=-1):
"""Three-state affine alignment. X consumes a; Y consumes b. Ties favor M,X,Y."""
if not a or not b or len(a) > 2000 or len(b) > 2000:
raise ValueError("Use two nonempty sequences, each at most 2,000 bases")
if gap_open > 0 or gap_extend > 0:
raise ValueError("Gap penalties must be nonpositive")
n, m = len(a), len(b)
scores = np.full((3,n+1,m+1), -np.inf)
trace = np.full((3,n+1,m+1), 255, dtype=np.uint8)
scores[0,0,0] = 0
if local:
scores[:] = 0
else:
for i in range(1,n+1):
scores[1,i,0] = gap_open + (i-1)*gap_extend
trace[1,i,0] = 0 if i == 1 else 1
for j in range(1,m+1):
scores[2,0,j] = gap_open + (j-1)*gap_extend
trace[2,0,j] = 0 if j == 1 else 2
best = (0.0,0,0,0)
for i in range(1,n+1):
for j in range(1,m+1):
values = scores[:,i-1,j-1]
parent = int(np.argmax(values))
scores[0,i,j] = values[parent] + (match if a[i-1] == b[j-1] else mismatch)
trace[0,i,j] = parent
values = (scores[0,i-1,j]+gap_open, scores[1,i-1,j]+gap_extend)
parent = int(np.argmax(values))
scores[1,i,j] = values[parent]; trace[1,i,j] = (0,1)[parent]
values = (scores[0,i,j-1]+gap_open, scores[2,i,j-1]+gap_extend)
parent = int(np.argmax(values))
scores[2,i,j] = values[parent]; trace[2,i,j] = (0,2)[parent]
if local:
for state in range(3):
if scores[state,i,j] <= 0:
scores[state,i,j] = 0; trace[state,i,j] = 255
state = int(np.argmax(scores[:,i,j]))
if scores[state,i,j] > best[0]:
best = (float(scores[state,i,j]),state,i,j)
if local:
score, state, i, j = best
else:
i, j = n, m; state = int(np.argmax(scores[:,i,j])); score = float(scores[state,i,j])
end = (i,j); aa, bb, path = [], [], [(i,j)]
while i or j:
if local and scores[state,i,j] == 0:
break
parent = int(trace[state,i,j])
if parent == 255:
raise ValueError("Invalid traceback")
if state == 0:
aa.append(a[i-1]); bb.append(b[j-1]); i-=1; j-=1
elif state == 1:
aa.append(a[i-1]); bb.append("-"); i-=1
else:
aa.append("-"); bb.append(b[j-1]); j-=1
state = parent; path.append((i,j))
return dict(a="".join(reversed(aa)),b="".join(reversed(bb)),score=score,
start=(i,j),end=end,path=list(reversed(path)),matrix=np.max(scores,axis=0))
print("Method ready: affine_align")
Method ready: affine_align
def alignment_summary(result, length_a, length_b):
a, b = result["a"], result["b"]
columns = len(a)
matches = sum(x == y and x != "-" for x,y in zip(a,b))
nongap = sum(x != "-" and y != "-" for x,y in zip(a,b))
return dict(score=result["score"],columns=columns,matches=matches,
identity_all_columns=matches/columns if columns else 0,
identity_nongap_columns=matches/nongap if nongap else 0,
coverage_a=(len(a)-a.count("-"))/length_a,
coverage_b=(len(b)-b.count("-"))/length_b,
gap_characters=a.count("-")+b.count("-"))
print("Method ready: alignment_summary")
Method ready: alignment_summary
Steps
1. Choose the sequence and scoring system
Align annotated coding sequences rather than mixing UTRs with CDS. Global alignment accounts for every base; local alignment selects a high-scoring segment. Match +2, mismatch −1, gap-open −3 and gap-extension −1 are transparent teaching parameters, not an inferred evolutionary model. A gap of length k costs open + (k−1)×extension. Both end gaps are penalized in our global model.
records = {r["gene"]:r for r in snapshot("sequences")["records"]}
sequences = {g:r["sequence"][r["cds_start"]-1:r["cds_end"]] for g,r in records.items()}
a,b = sequences["HBB"],sequences["HBD"]
SCORING = dict(match=2,mismatch=-1,gap_open=-3,gap_extend=-1)
print("HBB / HBD CDS lengths:",len(a),len(b)); print("Scoring:",SCORING)
HBB / HBD CDS lengths: 444 444
Scoring: {'match': 2, 'mismatch': -1, 'gap_open': -3, 'gap_extend': -1}
2. Run the three-state dynamic program
The function above keeps three states: M ends in a paired column, X ends in a gap in sequence b, and Y ends in a gap in sequence a. Extending a gap stays in the same gap state; opening one moves from M. Traceback follows stored parent states. Local alignment can restart at zero. Its coordinates refer to the original CDS, not the cropped output. Inspect the function before asking an AI assistant to modify it.
global_result = affine_align(a,b,**SCORING)
local_result = affine_align(a,b,local=True,**SCORING)
summary = pd.DataFrame([alignment_summary(r,len(a),len(b)) for r in (global_result,local_result)],
index=["global","local"])
display(summary)
print("Local 0-based intervals:",local_result["start"],local_result["end"])
| score | columns | matches | identity_all_columns | identity_nongap_columns | coverage_a | coverage_b | gap_characters | |
|---|---|---|---|---|---|---|---|---|
| global | 793.0 | 445 | 414 | 0.930337 | 0.934537 | 1.0 | 1.0 | 2 |
| local | 793.0 | 445 | 414 | 0.930337 | 0.934537 | 1.0 | 1.0 | 2 |
Local 0-based intervals: (0, 0) (444, 444)
3. Read a colored alignment
Green columns match; orange columns contain substitutions; gray indicates a gap. Print manageable 60-column blocks rather than one unreadable line. A visually similar alignment is not a gene annotation or a statistical homology test. Identity over all columns and identity over nongap paired columns answer different questions.
import html
blocks=[]
for start in range(0,len(global_result["a"]),60):
x,y=global_result["a"][start:start+60],global_result["b"][start:start+60]
colors=["#daf0e6" if p==q else "#eee" if "-" in (p,q) else "#ffe3cb" for p,q in zip(x,y)]
rows=["".join(f'<span style="background:{c};padding:1px">{html.escape(base)}</span>'
for base,c in zip(seq,colors)) for seq in (x,y)]
blocks.append(f'<p>Columns {start+1}–{start+len(x)}</p><div>HBB {rows[0]}</div><div>HBD {rows[1]}</div>')
display(HTML('<div style="font:12px/1.8 monospace;overflow-x:auto">'+"".join(blocks)+'</div>'))
Columns 1–60
Columns 61–120
Columns 121–180
Columns 181–240
Columns 241–300
Columns 301–360
Columns 361–420
Columns 421–445
4. Draw a sequence dot plot
Each dot marks an exact matching 9-mer. Coordinates are 0-based CDS starts. Repeated motifs can create off-diagonal points. This is a separate visualization of sequence similarity; it does not reuse the dynamic-programming traceback.
KMER=9
locations={}
for j in range(len(b)-KMER+1): locations.setdefault(b[j:j+KMER],[]).append(j)
dots=[(i,j) for i in range(len(a)-KMER+1) for j in locations.get(a[i:i+KMER],[])]
fig,ax=plt.subplots(figsize=(5,5))
if dots:
xy=np.array(dots); ax.scatter(xy[:,0],xy[:,1],s=8,color="#18785f",alpha=.6)
ax.set(xlabel="HBB CDS start (0-based)",ylabel="HBD CDS start (0-based)",title="Exact 9-mer matches")
plt.tight_layout();plt.show();print("Matching k-mer pairs:",len(dots))
Matching k-mer pairs: 268
5. Follow the traceback through the matrix
The heatmap is the maximum of M/X/Y at each cell. The blue path is a complete traceback from the global terminal state. It is not necessarily the only optimal path: ties follow a documented deterministic preference. Compare its shape with the dot plot, then change one scoring parameter at a time.
fig,ax=plt.subplots(figsize=(6,5))
matrix=global_result["matrix"].copy();matrix[~np.isfinite(matrix)]=np.nan
im=ax.imshow(matrix,origin="lower",aspect="auto",cmap="YlGn")
path=np.array(global_result["path"])
ax.plot(path[:,1],path[:,0],color="#3565af",lw=1.3)
ax.set(xlabel="HBD bases consumed",ylabel="HBB bases consumed",title="Global affine alignment and traceback")
fig.colorbar(im,ax=ax,label="Best state score");plt.tight_layout();plt.show()
6. Compare scoring sensitivity
Change the mismatch cost while keeping the input sequences and gap penalties fixed. Compare scores, gaps and all-column identity. These results answer a parameter-sensitivity question, not a hypothesis test for homology. Different scoring systems cannot be compared as though their raw scores used one common scale.
sensitivity=pd.DataFrame([dict(mismatch=penalty,**alignment_summary(affine_align(a,b,match=2,mismatch=penalty,gap_open=-3,gap_extend=-1),len(a),len(b))) for penalty in [-1,-2,-3]])
fig,axes=plt.subplots(1,2,figsize=(8,3))
axes[0].bar(sensitivity.mismatch.astype(str),sensitivity.identity_all_columns,color="#18785f")
axes[0].set(xlabel="Mismatch penalty",ylabel="Identity / all columns",ylim=(0,1),title="Identity under alternative costs")
axes[1].bar(sensitivity.mismatch.astype(str),sensitivity.gap_characters,color="#416aa6")
axes[1].set(xlabel="Mismatch penalty",ylabel="Gap characters",title="Gap sensitivity")
plt.tight_layout();plt.show();display(sensitivity[["mismatch","score","gap_characters"]])
| mismatch | score | gap_characters | |
|---|---|---|---|
| 0 | -1 | 793.0 | 2 |
| 1 | -2 | 764.0 | 2 |
| 2 | -3 | 738.0 | 6 |
7. Check coverage and export differences
Removing gaps must recover the original sequence for a global alignment. Differences are indexed by alignment column, not chromosome position. Gaps make those coordinates diverge. Keep the score parameters and accessions beside the output; an alignment without them cannot be reproduced reliably.
assert global_result["a"].replace("-","")==a
assert global_result["b"].replace("-","")==b
assert summary.loc["global","coverage_a"]==1 and summary.loc["global","coverage_b"]==1
differences=pd.DataFrame([dict(alignment_column=i+1,HBB=x,HBD=y)
for i,(x,y) in enumerate(zip(global_result["a"],global_result["b"])) if x!=y])
display(differences.head(12))
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
summary.to_csv(OUTPUT_DIR/"alignment-summary.csv")
differences.to_csv(OUTPUT_DIR/"alignment-differences.csv",index=False)
(OUTPUT_DIR/"alignment.txt").write_text(str(SCORING)+"\nHBB " +global_result["a"]+"\nHBD "+global_result["b"]+"\n")
print("Checks passed; differing columns:",len(differences))
| alignment_column | HBB | HBD | |
|---|---|---|---|
| 0 | 28 | T | A |
| 1 | 33 | C | T |
| 2 | 36 | T | C |
| 3 | 38 | C | A |
| 4 | 54 | G | A |
| 5 | 68 | A | C |
| 6 | 93 | G | A |
| 7 | 94 | C | T |
| 8 | 96 | G | A |
| 9 | 151 | A | T |
| 10 | 201 | A | G |
| 11 | 207 | C | A |
Checks passed; differing columns: 31
Exercises and reference explanation
- Calculate the penalty of a one-base and four-base gap with the default scoring.
- Rerun with gap-open −5, keeping every other parameter fixed. Compare scores and gap characters.
Reference explanation. The penalties are −3 and −6. A changed score can alter the chosen path. It does not mean the input sequences have changed. Removing global gaps must still recover each complete CDS.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Add an independent function that rescales an existing alignment by iterating over matches, mismatches and consecutive gap runs; compare its score to the dynamic-programming result.
Prompt 2
Plot global identity and coverage for three mismatch penalties. Use the existing affine implementation and keep all-column versus nongap identity distinct.
Human acceptance. Use synthetic short strings only as tests: identical sequences, a four-base deletion and a local shared segment. Validate score, recovered global strings and local source intervals.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NM_000518.5&rettype=fasta&retmode=text',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
Sources, snapshots and next steps
- RefSeq HBB/HBD — NCBI public sequence data; retrieved 2026-10-09. SHA-256
ea64e38e8e5ffad06f6d8770eaee066ef556e1488c4b9a6b340cbd445abe04ad. GenBank XML CDS uses 1-based inclusive coordinates
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.