从 FASTA 到基因与蛋白质

沿着真实 HBB 与 HBD 转录本提取编码区，将翻译结果与 RefSeq 注释逐字核对。

生物信息学 · 第 01 课🔋 离线核心分析⏱ 45 分钟内嵌真实数据📦 numpy · pandas · matplotlib · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python

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='01-fasta-to-protein'
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

读取内嵌快照

In [2] · Python

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

方法函数

In [3] · Python

def parse_fasta(text):
    records, label, chunks = [], None, []
    for line in text.splitlines():
        line = line.strip()
        if not line:
            continue
        if line.startswith(">"):
            if label is not None:
                records.append((label, "".join(chunks).upper()))
            label, chunks = line[1:].strip(), []
            if not label:
                raise ValueError("FASTA header is empty")
        else:
            if label is None:
                raise ValueError("FASTA requires a header before sequence lines")
            chunks.append(line)
    if label is not None:
        records.append((label, "".join(chunks).upper()))
    if not records:
        raise ValueError("No FASTA records")
    if len({label.split()[0] for label, seq in records}) != len(records):
        raise ValueError("Duplicate FASTA identifiers")
    for label, seq in records:
        if not seq or set(seq) - set("ACGTRYSWKMBDHVN"):
            raise ValueError("Empty or invalid DNA sequence: " + label)
    return records

print("Method ready: parse_fasta")

Method ready: parse_fasta

In [4] · Python

def translate_cds(sequence):
    sequence = sequence.upper()
    if len(sequence) % 3:
        raise ValueError("Complete CDS length must be divisible by three")
    bases = "TCAG"
    amino_acids = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
    codons = [a+b+c for a in bases for b in bases for c in bases]
    table = dict(zip(codons, amino_acids))
    protein = "".join(table.get(sequence[i:i+3], "X") for i in range(0,len(sequence),3))
    if "*" in protein[:-1]:
        raise ValueError("Internal stop codon: check the reading frame or genetic code")
    return protein.removesuffix("*")

print("Method ready: translate_cds")

Method ready: translate_cds

In [5] · Python

def gc_fraction(sequence):
    sequence = sequence.upper()
    denominator = sum(sequence.count(base) for base in "ACGT")
    if denominator == 0:
        return float("nan")
    return (sequence.count("G") + sequence.count("C")) / denominator

print("Method ready: gc_fraction")

Method ready: gc_fraction

分步分析

1. 先读记录，再读碱基

HBB 与 HBD 是人类旁系同源基因，不是两个患者的样本。RefSeq NM_000518.5 和 NM_000519.4 描述的是 mRNA，序列包含非翻译区（UTR）。两份记录的 CDS 均为 51..494。GenBank 采用从 1 开始、包含终点的坐标；Python 切片从 0 开始、不包含终点，因此应使用 [start-1:end]。保留 accession 的版本号，才能明确复现的是哪一条序列。

In [6] · Python

records = snapshot("sequences")["records"]
fasta = "\n".join(">"+r["accession"]+" "+r["gene"]+"\n"+r["sequence"] for r in records)
parsed = parse_fasta(fasta)
display(pd.DataFrame([dict(gene=r["gene"], accession=r["accession"],
    mrna_bases=len(r["sequence"]), cds=f'{r["cds_start"]}..{r["cds_end"]}',
    protein=r["protein_id"]) for r in records]))
print(fasta[:180] + "…")

	gene	accession	mrna_bases	cds	protein
0	HBB	NM_000518.5	628	51..494	NP_000509.1
1	HBD	NM_000519.4	620	51..494	NP_000510.1

>NM_000518.5 HBB
ACATTTGCTTCTGACACAACTGTGTTCACTAGCAACCTCAAACAGACACCATGGTGCATCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTGGTCTACCCTTG…

2. 绘制转录本结构

转录本示意图让我们在翻译前检查坐标。绿色矩形表示 CDS，灰色区域表示 UTR；横轴是 mRNA 坐标，不是染色体坐标。图形直接来自下载记录中的注释，不能据此推断基因的内含子结构。

In [7] · Python

fig, ax = plt.subplots(figsize=(8,3))
for i, r in enumerate(records):
    start, end, length = r["cds_start"], r["cds_end"], len(r["sequence"])
    ax.broken_barh([(1,start-1),(end+1,length-end)],(i-.22,.44),facecolors="#d9dedb")
    ax.broken_barh([(start,end-start+1)],(i-.22,.44),facecolors="#18785f")
    ax.text(start,i+.3,f"CDS {start}–{end}",fontsize=10)
ax.set(yticks=range(len(records)),yticklabels=[r["gene"] for r in records],
       xlabel="mRNA position (1-based)",title="Annotated coding regions and untranslated regions")
ax.spines[["top","right","left"]].set_visible(False)
plt.tight_layout(); plt.show()

3. 检查碱基组成与模糊碱基

GC 含量是描述序列组成的统计量。分母只计入 A/C/G/T，模糊碱基另行报告。两条转录本的 GC 差异不能直接证明表达量、稳定性或进化选择。后面的窗口参数可以修改，用来观察平滑程度如何改变图形。

In [8] · Python

from collections import Counter
composition = pd.DataFrame({r["gene"]:Counter(r["sequence"]) for r in records}).reindex(list("ACGT")).fillna(0)
composition.plot.bar(figsize=(7,3),color=["#18785f","#416aa6"],rot=0)
plt.ylabel("Base count"); plt.title("Composition of the complete mRNA sequences")
plt.tight_layout(); plt.show()
display(pd.DataFrame([dict(gene=r["gene"],gc_fraction=gc_fraction(r["sequence"]),
    ambiguous_bases=sum(b not in "ACGT" for b in r["sequence"])) for r in records]))

	gene	gc_fraction	ambiguous_bases
0	HBB	0.512739	0
1	HBD	0.506452	0

4. 比较局部 GC 含量

每个点使用 80 个碱基的窗口，步长为 20 个碱基，位置表示窗口中点；末端不足一个窗口的片段被省略。短转录本只支持少量窗口，不应将曲线视作全基因组规律。

In [9] · Python

WINDOW, STEP = 80, 20
fig, ax = plt.subplots(figsize=(8,3))
for r in records:
    starts = range(0,len(r["sequence"])-WINDOW+1,STEP)
    ax.plot([s+(WINDOW+1)/2 for s in starts],
            [100*gc_fraction(r["sequence"][s:s+WINDOW]) for s in starts],marker="o",label=r["gene"])
ax.set(xlabel="Window midpoint (mRNA coordinate)",ylabel="GC (%)",title=f"Local GC: {WINDOW}-base windows")
ax.legend(); plt.tight_layout(); plt.show()

5. 翻译注释中的编码序列

使用标准遗传密码表，不能简单地从整条转录本的第一个 ATG 开始翻译。完整 CDS 包含末端终止密码子，RefSeq 蛋白序列不包含该终止符。本例 444 个编码碱基对应 147 个氨基酸及一个终止密码子。翻译产物包含起始甲硫氨酸，而成熟血红蛋白链的编号可能不包含它。

In [10] · Python

proteins, cds_records = {}, {}
for r in records:
    cds = r["sequence"][r["cds_start"]-1:r["cds_end"]]
    protein = translate_cds(cds)
    assert protein == r["translation"], "Translation differs from RefSeq"
    proteins[r["gene"]], cds_records[r["gene"]] = protein, cds
    print(r["gene"], "CDS:",len(cds), "bases; protein:",len(protein),"aa; RefSeq match: True")
    print(protein)

HBB CDS: 444 bases; protein: 147 aa; RefSeq match: True
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH
HBD CDS: 444 bases; protein: 147 aa; RefSeq match: True
MVHLTPEEKTAVNALWGKVNVDAVGGEALGRLLVVYPWTQRFFESFGDLSSPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFSQLSELHCDKLHVDPENFRLLGNVLVCVLARNFGKEFTPQMQAAYQKVVAGVANALAHKYH

6. 校验并导出

断言可以使坐标或翻译错误立即显现。分别导出核酸和蛋白序列，方便在其他工具中继续分析。输出目录创建在 notebook 所在文件夹，重复运行只会更新本课生成的结果文件。

In [11] · Python

assert len(parsed) == 2
assert all(len(cds) == 444 for cds in cds_records.values())
assert all(len(p) == 147 for p in proteins.values())
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
(OUTPUT_DIR/"transcripts.fasta").write_text(fasta+"\n")
(OUTPUT_DIR/"proteins.fasta").write_text("".join(f">{g}\n{p}\n" for g,p in proteins.items()))
pd.DataFrame([dict(gene=g,cds_bases=len(cds_records[g]),protein_aa=len(p),gc=gc_fraction(cds_records[g]))
              for g,p in proteins.items()]).to_csv(OUTPUT_DIR/"sequence-summary.csv",index=False)
print("Checks passed. Exported:",[p.name for p in OUTPUT_DIR.iterdir()])

Checks passed. Exported: ['proteins.fasta', 'transcripts.fasta', 'sequence-summary.csv']

自主练习与参考答案

解释为什么 CDS 51..494 应使用 [50:494]，并预测把起始索引改成 51 会怎样。
将 WINDOW 从 80 改为 40，保持原始转录本不变，更新图中的窗口说明。

参考解释。 正确 CDS 为 444 个碱基。使用起始索引 51 会删去一个碱基，长度变成 443，无法通过完整 CDS 校验。较小窗口展示更多局部变化，不能据此证明调控机制。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

只使用标准库及当前包，在每个 GC 值旁添加模糊碱基计数，明确保留 A/C/G/T 分母。

提示词 2

添加支持 IUPAC DNA 符号的反向互补函数，并断言执行两次后恢复原序列。解释为何对注释中的正向 CDS 做反向互补不能代替链方向注释。

人工验收。 检查两份 RefSeq 翻译仍完全一致；N-only GC 应为 NaN，非法符号应报错，两次反向互补应恢复原序列。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [12] · Python

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

来源、快照与下一步

RefSeq HBB/HBD — NCBI public sequence data; retrieved 2026-10-09. SHA-256 ea64e38e8e5ffad06f6d8770eaee066ef556e1488c4b9a6b340cbd445abe04ad. GenBank XML CDS uses 1-based inclusive coordinates

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。