比对 β 与 δ 珠蛋白基因

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

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']

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

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}

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"])

Local 0-based intervals: (0, 0) (444, 444)

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>'))

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

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()

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

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))

Checks passed; differing columns: 31

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