实验 1：序列基础（FASTA · GC 含量 · 反向互补 · 翻译 · ORF）

预计用时：6–10 小时，建议分 3 次做完。
对应学习指南「项目一：DNA 序列基础分析」，另外加上了翻译和 ORF，用真实的人类 TP53 mRNA 收尾。
开始前先读同目录的 前置知识.md。

一、目标

做完这个实验，你应该能：

不靠库，手写一个 FASTA 解析器，并说清楚 FASTA 格式有哪些坑
区分格式问题（小写、空格、换行）和内容问题（非法字符、DNA/RNA 混用），分开处理
正确计算 GC 含量：知道分母为什么只能算确定碱基
对 DNA、RNA、简并碱基都能正确做反向互补
用密码表翻译序列，在两条链 × 三个读框里找 ORF
从真实的 TP53 mRNA 里找出 p53 蛋白的编码区，和 UniProt 的数据对上

二、前置

类别	需要什么	在哪补
生物	中心法则、碱基配对、反向平行、密码子、ORF	`前置知识.md` 第一部分
Python	字符串方法、dict / Counter、列表推导、dataclass、异常	`前置知识.md` 第二部分
工具	会用 `course start / check / run`	课程根目录 `README.md`

三、数据

文件	内容	来源
`数据/训练序列.fasta`	6 条短序列，故意混进小写、空格、空行、RNA、N、X、简并码	教学构造
`数据/TP53_NM_000546.6.fasta`	人 TP53 基因 mRNA 全长，共 2512 nt	NCBI RefSeq `NM_000546.6`（2026-09-27 通过 Entrez efetch 下载）

想自己去 NCBI 取序列，可以做选做题 1。

四、任务

先运行 course start 1，然后在 交付/seqtools.py 里按顺序完成。每个任务都有对应的测试，测试名以 任务N_ 开头，所以可以用 course check 1 -k 任务3 只测一个任务。

任务 1：手写 FASTA 解析器 `read_fasta`

要求：标题行拆成 id 和 description；多行序列要拼接；空行跳过；标题行之前就出现序列时报错。
为什么不直接用 Bio.SeqIO：以后遇到格式有点问题的文件（这种文件很常见），你得知道解析器在哪一步做了什么假设。测试里会拿 SeqIO 的结果来对你的答案。

任务 2：清洗与校验 `clean` / `molecule_type` / `invalid_chars`

清洗只管格式：转大写、去掉非字母。X 这类非法字母要保留下来。
校验只管内容：按 IUPAC 字母表找出非法字符。N 表示"不确定是哪个碱基"，是合法的；X 不是核酸字符。
T 和 U 同时出现时直接报错：这说明数据本身有问题，程序不应该自己猜是哪一种。

任务 3：组成统计 `base_counts` / `gc_content` / `gc_windows`

GC% = (G+C) / 确定碱基数 × 100。N、简并码、非法字符都不进分母。
一个确定碱基都没有时返回 nan，不要返回 0：0 表示"GC 含量是 0"，而 nan 表示"没法算"，两者意思不同。
滑窗 GC：后面画 TP53 的 GC 分布图要用。

任务 4：反向互补 `reverse_complement`

DNA 和 RNA 用各自的互补表：RNA 里 A 和 U 配对，结果不能出现 T。
简并码也有互补，例如 R（A 或 G）的互补是 Y（C 或 T）。
碰到没法互补的字符，要抛出 ValueError，不要偷偷换成 N。一个看起来正常的错误结果，比直接报错更难发现。

任务 5：翻译与 ORF `translate` / `find_orfs`

翻译：U 先换成 T，再查 CODON_TABLE。含 N 的密码子翻成 X，终止密码子翻成 *。
ORF：两条链 × 三个读框，从 ATG 开始，读到同一读框的第一个终止密码子。
坐标是最容易错的地方。负链上找到的 ORF 要换算回正链坐标，这一点在函数的 docstring 里写明了。
验收：TP53 mRNA 里最长的 ORF 应该是 393 aa，以 MEEPQSDPSV 开头。这就是 p53 蛋白（UniProt P04637）。

任务 6：汇总与交付 `summarize`，然后 `course run 1`

汇总表每条序列一行。含非法字符的序列，反向互补那一列留空。
运行 course run 1 后，交付/ 里会出现三个文件：序列汇总.csv、TP53_最长ORF.txt、TP53_GC滑窗.png。

思考题（写进 `交付/思考题.md`，每题 2–4 句）

sample4 按旧算法得到 GC% = 30.3%，按新算法是 50%。哪个对？在真实的基因组分析里，这种误差会造成什么后果？
用 gc_content 分别算 5'UTR（ORF 之前）、编码区、3'UTR（ORF 之后）三段的 GC%，再对照 TP53_GC滑窗.png：哪一段最高，哪一段最低？查资料解释可能的原因（关键词：5'UTR 二级结构、3'UTR 富含 AU 的元件）。
find_orfs 在 TP53 mRNA 上还找到了其他比较短的 ORF，它们是真的基因吗？只靠 ORF 预测基因，有什么局限？

选做

用 Bio.Entrez.efetch 自己下载 NM_000546.6，确认和 数据/ 里的文件一模一样。记得填 Entrez.email。
在 Rosalind 上做 DNA、RNA、REVC、GC、PROT、ORF 这 6 题，正好对应任务 2–5。

五、验收标准

检查点	通过标准	对应测试
解析	6 条记录；多行拼接、空行跳过；结果和 Biopython 一致	`任务1_*`
校验	sample4 的非法字符只有 `X`；T 和 U 混用时报错	`任务2_*`
GC	sample1 = 47.5；sample4 = 50.0；全 N 返回 nan	`任务3_*`
反向互补	sample1 = `ACGTACGATGACTAGCATACGCATGCATGCATGCATGCAT`；RNA 结果里没有 T；和 Biopython 一致	`任务4_*`
翻译 / ORF	sample3 = `MAMAGMPCMAMPC`；TP53 最长 ORF 为 393 aa，坐标 143..1324	`任务5_*`
交付	`course check 1` 全部通过，`交付/` 里有 3 个文件和思考题	`任务6_*`

六、常见错误

现象	原因	解决
sample4 的 GC 算出 30.3	分母把 N 和 X 也算进去了	分母只数 A/C/G/T/U
RNA 的反向互补结果里出现 T	用了 DNA 的互补表	先用 `molecule_type` 判断类型，再选互补表
`clean` 之后中文描述混进了序列	中文字符的 `isalpha()` 也返回 True	再加一个 `isascii()` 判断
负链 ORF 坐标对不上	直接用了反向互补链上的下标	正链坐标 = `len - e` 到 `len - s`
TP53 找到的最长 ORF 不是 393	同一读框里遇到新的 ATG 就重新开始了	已经在 ORF 里时，忽略后面的 ATG
翻译结果比预期多一个 X	末尾不足 3 个碱基的部分也翻译了	`range(0, len - len % 3, 3)`
测试报 `NotImplementedError`	这个函数还没写	看测试输出最后的"还没写的 TODO"清单
读文件报 `UnicodeDecodeError`	没有指定编码	`open(path, encoding="utf-8")`

七、和旧版的区别（为什么重做）

旧版在 归档/v1_实验1_序列分析_表达可视化/，有 4 个实打实的错误：

GC 分母错了：X 和 N 被算进分母，sample4 算成 30.3%，正确值是 50%。
RNA 反向互补错了：sample2 的结果里出现了 T。
非法字符没有处理：说明里写了"过滤"，代码里却原样保留，还把含 X 的序列拿去反向互补。
模板本身有语法错误：raw.upper()'''???''' 和 df.???(...)，根本跑不起来。

旧版的表达可视化部分已经换成了实验 2：真实数据加上正规的差异表达分析。

实验 1 前置知识：生物概念 ↔ 编程实现对照表

给没有生物背景的同学看。每个生物概念都配了它在代码里的样子。
做实验时遇到不懂的概念，回来查这张表就行，不用先背下来。

一、生物概念 ↔ 代码

生物概念	一句话解释	在代码里	本实验哪里用到
DNA	由 A/C/G/T 四种碱基组成的双链分子，存储遗传信息	只含 `ACGT` 的字符串	全部
RNA	DNA 转录出的单链拷贝，T 换成了 U	只含 `ACGU` 的字符串	`molecule_type`、`reverse_complement`
碱基配对	A 配 T（RNA 里 A 配 U），G 配 C	一个查找表：dict `{"A": "T", ...}`	任务 4
反向平行	双链的两条链方向相反，一条是 5'→3'，另一条是 3'→5'	所以求另一条链要"互补 + 倒序"：`seq[::-1]`	任务 4
反向互补	把另一条链按 5'→3' 方向写出来	`"".join(table[b] for b in reversed(seq))`	任务 4、5
IUPAC 简并码	测序时不能确定是哪个碱基，就用一个字母表示几种可能：N=任意，R=A/G，Y=C/T…	合法字符，但不算"确定碱基"	任务 2、3
GC 含量	G 和 C 占碱基的比例。GC 配对有 3 个氢键，所以 GC 高的片段更稳定，PCR 引物设计要参考它	`(G+C) / 确定碱基数 × 100`	任务 3
密码子	3 个碱基编码 1 个氨基酸。共 64 种，其中 3 种是"终止"	按 3 个一组切片 + 查字典 `CODON_TABLE`	任务 5
读框（frame）	从第 0、1、2 位开始切三联体，读出来的蛋白完全不同	`range(frame, n - 2, 3)`	任务 5
起始 / 终止密码子	ATG 开始翻译（对应甲硫氨酸 M）；TAA、TAG、TGA 结束翻译	`START_CODON`、`STOP_CODONS`	任务 5
ORF（开放读框）	从 ATG 到同读框第一个终止密码子之间的一段，可能编码蛋白	一个状态机："已经在 ORF 里 / 不在 ORF 里"	任务 5
mRNA 的结构	5'UTR（不翻译）+ CDS（编码区）+ 3'UTR（不翻译）	三个切片 `seq[:s]`、`seq[s:e]`、`seq[e:]`	任务 5、思考题
RefSeq 编号	NCBI 给参考序列的 ID。NM_ 开头是 mRNA，`.6` 是版本号	FASTA 标题行里的 id	数据

二、Python 前置

知识点	最小例子	用在
逐行读文件	`with open(p, encoding="utf-8") as f: for line in f:`	任务 1
字符串方法	`s.upper()`、`s.strip()`、`s.startswith(">")`、`s.split(maxsplit=1)`	任务 1、2
字符判断	`"x".isalpha()`、`"中".isascii()` 返回 False	任务 2
集合运算	`set(seq) - alphabet` 得到"不在字母表里的字符"	任务 2
计数	`Counter("AACG")` 的结果是 `{"A": 2, "C": 1, "G": 1}`	任务 3
抛异常	`raise ValueError("说清楚哪里错了")`	任务 1、2、4
切片与步长	`s[i:i+3]`、`range(0, len(s), 3)`、`reversed(s)`	任务 4、5
dataclass	`ORF("+", 0, 2, 38, "MAAA")`，用 `o.start` 取字段	任务 5
pandas 建表	`pd.DataFrame([{"id": ..., "gc": ...}, ...])`	任务 6
nan	`math.nan`；判断要用 `math.isnan(x)`，因为 `nan == nan` 是 False	任务 3

三、脏数据清单（训练序列里都埋了）

情况	例子	属于	怎么处理
小写	`atgc`	格式	`clean` 转大写
空格、制表符、数字	`ATG CGT`、`12 ATG`	格式	`clean` 去掉
空行、跨行	一条序列分成几行写	格式	`read_fasta` 拼接、跳过
N	`ATGNNATG`	内容，合法	保留，不进 GC 的分母
简并码	`ACGTRY`	内容，合法	保留，反向互补时用简并码互补表
非核酸字母	`X`	内容，非法	`invalid_chars` 报出来；不做反向互补
T 和 U 混用	`ACGTU`	内容，非法	`molecule_type` 报错

参考实现：seqtools

卡住超过 20 分钟再看。对照着看自己是哪一步想岔了 —— 抄一遍没用。

展开参考解

看答案

"""实验 1 · 序列基础（参考解）

先自己在 交付/seqtools.py 里写，卡住超过 20 分钟再来对照。
每个函数的"为什么这样做"写在 实验说明.md 对应的任务里。
"""
from __future__ import annotations

import math
from collections import Counter
from dataclasses import dataclass
from pathlib import Path

from Bio.Data import CodonTable

HERE = Path(__file__).resolve().parent
DATA = HERE.parent / "数据"

# IUPAC 核酸字母表：4 个确定碱基 + 11 个简并码（N = 任意碱基）
DEFINITE_DNA = set("ACGT")
DEFINITE_RNA = set("ACGU")
AMBIGUOUS = set("RYSWKMBDHVN")
IUPAC_DNA = DEFINITE_DNA | AMBIGUOUS
IUPAC_RNA = DEFINITE_RNA | AMBIGUOUS

# 互补：确定碱基两两配对；简并码也有互补（R=A/G 的互补是 Y=C/T，依此类推）
_AMBIG_COMPLEMENT = {"R": "Y", "Y": "R", "S": "S", "W": "W", "K": "M", "M": "K",
                     "B": "V", "V": "B", "D": "H", "H": "D", "N": "N"}
DNA_COMPLEMENT = {"A": "T", "T": "A", "G": "C", "C": "G", **_AMBIG_COMPLEMENT}
RNA_COMPLEMENT = {"A": "U", "U": "A", "G": "C", "C": "G", **_AMBIG_COMPLEMENT}

# 标准遗传密码表（NCBI 1 号表），DNA 写法。终止密码子翻成 "*"。
_STANDARD = CodonTable.unambiguous_dna_by_id[1]
CODON_TABLE: dict[str, str] = {**_STANDARD.forward_table, **{c: "*" for c in _STANDARD.stop_codons}}
START_CODON = "ATG"
STOP_CODONS = set(_STANDARD.stop_codons)


@dataclass
class FastaRecord:
    id: str
    description: str
    seq: str  # 原始序列：多行已拼接，但还没清洗（可能有小写、空格）


@dataclass
class ORF:
    strand: str   # "+" 正链，"-" 反向互补链
    frame: int    # 0/1/2：从该链第几个碱基开始按三联体读
    start: int    # 在【原始正链】上的 0 起始坐标，左闭
    end: int      # 在【原始正链】上的坐标，右开（包含终止密码子）
    protein: str  # 翻译出的蛋白（不含终止符 *）


# ---------------------------------------------------------------- 任务 1：读 FASTA
def read_fasta(path) -> list[FastaRecord]:
    records: list[FastaRecord] = []
    header = None
    chunks: list[str] = []
    with open(path, encoding="utf-8") as f:
        for lineno, line in enumerate(f, 1):
            line = line.rstrip("\r\n")
            if not line.strip():
                continue
            if line.startswith(">"):
                if header is not None:
                    records.append(_make_record(header, chunks))
                header, chunks = line[1:].strip(), []
            else:
                if header is None:
                    raise ValueError(f"第 {lineno} 行：在第一个 '>' 标题行之前出现了序列")
                chunks.append(line.strip())
    if header is not None:
        records.append(_make_record(header, chunks))
    return records


def _make_record(header: str, chunks: list[str]) -> FastaRecord:
    parts = header.split(maxsplit=1)
    return FastaRecord(id=parts[0], description=parts[1] if len(parts) > 1 else "", seq="".join(chunks))


# ---------------------------------------------------------------- 任务 2：清洗与校验
def clean(raw: str) -> str:
    """转大写，去掉空白、数字等所有非字母字符。非法字母（如 X）保留，交给 invalid_chars 报告。"""
    return "".join(ch for ch in raw.upper() if ch.isascii() and ch.isalpha())


def molecule_type(seq: str) -> str:
    has_t, has_u = "T" in seq, "U" in seq
    if has_t and has_u:
        raise ValueError("序列同时含 T 和 U，不是合法的 DNA 或 RNA")
    return "RNA" if has_u else "DNA"


def invalid_chars(seq: str, molecule: str | None = None) -> list[str]:
    molecule = molecule or molecule_type(seq)
    alphabet = IUPAC_RNA if molecule == "RNA" else IUPAC_DNA
    return sorted(set(seq) - alphabet)


# ---------------------------------------------------------------- 任务 3：组成统计
def base_counts(seq: str) -> dict[str, int]:
    return dict(Counter(seq))


def gc_content(seq: str) -> float:
    """GC% = (G+C) / 确定碱基数 × 100。N、简并码、非法字符都不进分母；没有确定碱基返回 nan。"""
    counts = Counter(seq)
    definite = sum(counts[b] for b in "ACGTU")
    if definite == 0:
        return math.nan
    return (counts["G"] + counts["C"]) / definite * 100


def gc_windows(seq: str, window: int = 100, step: int = 25) -> list[tuple[int, float]]:
    if window <= 0 or step <= 0:
        raise ValueError("window 和 step 必须为正")
    return [(i, gc_content(seq[i:i + window])) for i in range(0, len(seq) - window + 1, step)]


# ---------------------------------------------------------------- 任务 4：反向互补
def reverse_complement(seq: str, molecule: str | None = None) -> str:
    molecule = molecule or molecule_type(seq)
    table = RNA_COMPLEMENT if molecule == "RNA" else DNA_COMPLEMENT
    bad = sorted(set(seq) - set(table))
    if bad:
        raise ValueError(f"无法互补的字符：{''.join(bad)}（先清洗/过滤再做反向互补）")
    return "".join(table[b] for b in reversed(seq))


# ---------------------------------------------------------------- 任务 5：翻译与 ORF
def translate(seq: str, to_stop: bool = False) -> str:
    dna = seq.replace("U", "T")
    protein = []
    for i in range(0, len(dna) - len(dna) % 3, 3):
        aa = CODON_TABLE.get(dna[i:i + 3], "X")  # 含 N 等简并码的密码子无法确定 → X
        if to_stop and aa == "*":
            break
        protein.append(aa)
    return "".join(protein)


def find_orfs(seq: str, min_protein_len: int = 30) -> list[ORF]:
    """两条链 × 三个读框里找 ATG…终止 的开放读框；同一读框里取最靠前的 ATG（最长）。"""
    dna = seq.replace("U", "T")
    n = len(dna)
    rc = reverse_complement(dna, "DNA")
    orfs: list[ORF] = []
    for strand, s in (("+", dna), ("-", rc)):
        for frame in range(3):
            start = None
            for i in range(frame, n - 2, 3):
                codon = s[i:i + 3]
                if start is None and codon == START_CODON:
                    start = i
                elif start is not None and codon in STOP_CODONS:
                    protein = translate(s[start:i])
                    if len(protein) >= min_protein_len:
                        end = i + 3
                        fwd = (start, end) if strand == "+" else (n - end, n - start)
                        orfs.append(ORF(strand, frame, fwd[0], fwd[1], protein))
                    start = None
    orfs.sort(key=lambda o: (-len(o.protein), o.start))
    return orfs


# ---------------------------------------------------------------- 任务 6：汇总与交付
def summarize(records: list[FastaRecord]):
    import pandas as pd

    rows = []
    for r in records:
        seq = clean(r.seq)
        if not seq:
            continue
        mol = molecule_type(seq)
        bad = invalid_chars(seq, mol)
        counts = base_counts(seq)
        rows.append({
            "id": r.id,
            "molecule": mol,
            "length": len(seq),
            **{b: counts.get(b, 0) for b in "ACGTU"},
            "N": counts.get("N", 0),
            "ambiguous": sum(v for k, v in counts.items() if k in AMBIGUOUS - {"N"}),
            "invalid": "".join(bad),
            "gc_percent": round(gc_content(seq), 2),
            # 有非法字符时不给反向互补：宁可空着，也不输出一条看起来像真的错误序列
            "reverse_complement": "" if bad else reverse_complement(seq, mol),
        })
    return pd.DataFrame(rows)


def main(out_dir) -> None:
    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)

    table = summarize(read_fasta(DATA / "训练序列.fasta"))
    table.to_csv(out / "序列汇总.csv", index=False, encoding="utf-8-sig")
    print(table.drop(columns="reverse_complement").to_string(index=False))

    tp53 = read_fasta(DATA / "TP53_NM_000546.6.fasta")[0]
    seq = clean(tp53.seq)
    best = find_orfs(seq, min_protein_len=100)[0]
    (out / "TP53_最长ORF.txt").write_text(
        f"来源：{tp53.id} {tp53.description}\n"
        f"mRNA 长度 {len(seq)} nt，全长 GC% {gc_content(seq):.1f}\n"
        f"最长 ORF：{best.strand} 链，读框 {best.frame}，坐标 {best.start + 1}..{best.end}（1 起始，含终止密码子）\n"
        f"蛋白长度 {len(best.protein)} aa\n\n{best.protein}\n",
        encoding="utf-8",
    )
    print(f"\nTP53 最长 ORF：{best.start + 1}..{best.end}，{len(best.protein)} aa，开头 {best.protein[:10]}")

    windows = gc_windows(seq, window=100, step=10)
    xs = [i + 50 for i, _ in windows]
    ys = [g for _, g in windows]
    fig, ax = plt.subplots(figsize=(9, 3.5))
    ax.plot(xs, ys, lw=1.2)
    ax.axvspan(best.start, best.end, color="orange", alpha=0.15, label="longest ORF (CDS)")
    ax.axhline(gc_content(seq), ls="--", color="gray", lw=0.8, label="whole-mRNA GC%")
    ax.set(xlabel="position (nt)", ylabel="GC% (100-nt window)", title=f"{tp53.id} TP53 mRNA")
    ax.legend(loc="lower left", fontsize=8)
    fig.tight_layout()
    fig.savefig(out / "TP53_GC滑窗.png", dpi=150)
    plt.close(fig)

目录实验 2：批量 RNA-seq 差异表达（真实数据 · PyDESeq2） →