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

预计用时：10–14 小时，建议分 4 次做完。
对应学习指南「项目二：基因表达数据可视化」。这次不再对着表达量排序看热图，而是回答一个真实的研究问题，
用的是 DESeq2 官方教程采用的经典数据集。
开始前先读同目录的 前置知识.md。

一、研究问题与数据

问题：糖皮质激素（地塞米松）是治疗哮喘的常用药。它作用于气道平滑肌细胞时，改变了哪些基因的表达？

数据：airway 数据集，出自 Himes et al., PLoS One 2014（GEO 编号 GSE52778）。

4 位供者的气道平滑肌细胞系，每个细胞系各有一份未处理、一份地塞米松 1 µM 处理 18 小时，共 8 个样本
计数来自 recount3：它把 SRA 上的公开数据用同一套流程重新比对过（STAR + GENCODE v26）
下载：course fetch 2，约 4 MB，下载后自动做 sha256 校验

这是"配对设计"：同一个细胞系，处理前后各测一次。这个设计决定了后面的统计模型怎么写。

二、目标

读懂真实数据：样本元数据、基因注释（GTF）、计数矩阵，并把它们对齐
理解为什么必须先过滤和标准化，再比较表达量
用 PCA 做质控，从图上看出"处理效应"和"细胞系效应"
用 PyDESeq2 做配对设计的差异表达分析，读懂 log2FC、p 值、padj
画火山图、MA 图、热图，并用已知的激素应答基因检验结果靠不靠谱
说清楚旧版实验的"均值 Top10 + 2 倍"错在哪里

三、前置

类别	需要什么	在哪补
生物	转录组、RNA-seq 测的是什么、基因 ID 与基因名、糖皮质激素受体	`前置知识.md` 第一部分
统计	为什么要标准化、p 值与多重检验、FDR、负二项分布（知道名字即可）	`前置知识.md` 第二部分
Python	pandas 的索引对齐、groupby、布尔筛选；numpy 的 SVD	`前置知识.md` 第三部分

四、任务

先运行 course fetch 2 和 course start 2，然后在 交付/de.py 里完成。DESeq2 相关的测试要跑 30 秒左右，前面几个任务可以单独测：course check 2 -k "任务1 or 任务2"。

任务 1：样本表 `parse_attributes` / `load_samples`

recount3 把细胞系和处理方式塞在一个字符串里：cell line;;N61311|...|treatment;;Untreated，需要拆开。
只保留"未处理"和"地塞米松"两组（原研究还有沙丁胺醇组，本实验不用）。
验收：8 个样本，4 个细胞系，每个细胞系正好一对。

任务 2：计数矩阵与注释 `coverage_to_counts` / `load_gene_annotation`

recount3 存的是碱基覆盖度之和，不是读数。除以平均比对长度，才约等于读数。这是 recount3 官方 compute_read_counts() 的做法。
GTF 是基因组注释的标准格式。ID 要去掉版本号，Y 染色体 PAR 区的重复拷贝要删掉。
验收：矩阵是 63811 × 8。FKBP5 在同一细胞系里的计数，未处理是 245，处理后是 4355。

任务 3：过滤与标准化 `filter_low_counts` / `cpm` / `log_cpm`

各样本测到的读数总量从 1900 万到 3200 万不等。直接比较原始计数，等于拿不同大小的碗比米粒数。
低表达基因的计数几乎全是噪声，还会拖累多重检验的校正，所以要先过滤掉。
验收：过滤后剩 19619 个基因；CPM 每一列加起来是 100 万。

任务 4：PCA `pca`

手写 PCA：中心化 → SVD → 得分。测试里会拿 sklearn 的结果对答案。
看图 course run 2 生成的 pca.png：PC1 把处理组和未处理组分开，PC2 大体按细胞系分开。这说明细胞系本身就是一个很大的差异来源，所以必须用配对设计把它扣掉。

任务 5：差异表达 `run_deseq` / `classify` / `naive_fold_change`

设计公式 ~cell + dex：先扣除细胞系的效应，再估计地塞米松的效应。
contrast 的方向是 dex 相对 untreated，写反了 FKBP5 就会变成"下调"。
显著的标准：padj < 0.05 且 |log2FC| ≥ 1。padj 已经对约两万次检验做过 BH 校正。
验收：6 个已知的激素应答基因（FKBP5、TSC22D3、PER1、ZBTB16、KLF15、DUSP1）全部显著上调。上调约 553 个、下调约 485 个，允许 ±5% 的误差。

任务 6：出图与报告：`course run 2`

交付/ 里会生成：pca.png、volcano.png、ma.png、heatmap_top30.png、de_results.csv、结果摘要.md。

思考题（写进 `交付/思考题.md`）

把设计公式改成 ~dex，也就是不考虑细胞系，重新跑一遍。显著基因变多了还是变少了？结合 PCA 图解释原因。
结果摘要.md 里写着：旧方法挑出了 825 个"两倍以上"的基因，其中 191 个在 DESeq2 里并不显著。从 de_results.csv 里挑两三个这样的基因，看看它们的 baseMean 和各样本的计数，说明旧方法为什么会上当。
为什么要看 padj，而不是 pvalue？如果对两万个基因都用 p < 0.05 作判断，大约会有多少个假阳性？
查一下 FKBP5 和 TSC22D3（也叫 GILZ）的功能，用一两句话说明它们为什么会被糖皮质激素上调。

选做

用 pydeseq2 的 DeseqStats.lfc_shrink() 做 log2FC 收缩，比较收缩前后的 MA 图。
把沙丁胺醇组也加进来，设计成 ~cell + treatment，比较两种药物的效应。

五、验收标准

检查点	通过标准	对应测试
样本表	8 个样本，4 个细胞系 × 2 种处理	`任务1_*`
计数	63811 × 8；SRR1039508 的总读数为 20913038；FKBP5 为 245 → 4355	`任务2_*`
过滤 / 标准化	剩 19619 个基因；CPM 每列之和为 1e6	`任务3_*`
PCA	和 sklearn 一致（只差正负号）；PC1 分开两组，约 40%	`任务4_*`
DESeq2	方向正确；6 个已知基因 padj < 1e-10 且 log2FC > 2；上调约 553、下调约 485	`任务5_*`
交付	`course check 2` 全部通过；6 个交付文件齐全；思考题已写	`任务6_*`

六、常见错误

现象	原因	解决
`course check 2` 全部跳过，提示缺数据	还没下载	`course fetch 2`
计数矩阵全是很大的数（上亿）	没把覆盖度换算成读数	除以 `avg_len`
注释 join 后大量 gene_name 是 NaN	一边的 ID 带版本号，另一边没有	两边都去掉 `.xx`
`index is not unique`	`_PAR_Y` 的重复没去掉	过滤掉以 `_PAR_Y` 结尾的 ID
PyDESeq2 报 counts 形状不对	传进去的是基因×样本	转置 `counts.T`
FKBP5 的 log2FC 是负数	contrast 方向写反了	`["dex", "dex", "untreated"]`
PCA 得分和 sklearn 差一个负号	主成分的方向本来就不唯一	正常，测试比较的是绝对值
运行时出现一条 "residual degrees of freedom" 警告	8 个样本、5 个参数，自由度只剩 3	已知情况，可以屏蔽（见模板说明）

七、这次和旧版的区别

旧版用的是 20 个基因 × 6 个样本的编造数据，做法有三个问题：

按"肿瘤组均值 Top10"挑基因，挑出来的只是表达量高的基因，不是差异表达的基因
用"> 2 倍"代替显著性检验：没有做标准化，没有检验，也没有做多重检验校正
结论写成"显著上调"，但数据里根本没有支持"显著"的统计证据

这一版用真实数据、正规流程，还有能拿文献核对的已知答案。

实验 2 前置知识：生物概念 ↔ 统计 ↔ 代码对照表

一、生物概念

概念	一句话解释	在数据/代码里
转录组	某一时刻细胞里所有 RNA 的集合，反映"哪些基因正在工作"	计数矩阵的每一列
RNA-seq	把 RNA 打碎、测序，再把读段（read）比对回基因组，数每个基因上落了多少读段	计数矩阵里的整数
读段计数（count）	落在某个基因上的读段数。基因越长、测序越深，计数越大	`counts.loc[基因, 样本]`
覆盖度（coverage）	每个碱基被多少读段覆盖。recount3 存的是它在整个基因上的总和	需要除以读长，换算成计数
Ensembl 基因 ID	`ENSG00000096060.14`：稳定的机器编号，`.14` 是版本号	计数矩阵的行索引
基因名（symbol）	`FKBP5`：人读的名字，可能改名，也可能重名	从 GTF 注释里查
GTF	基因组注释格式：每行是一个基因、转录本或外显子，第 9 列写属性	`load_gene_annotation`
糖皮质激素受体（GR）	地塞米松进入细胞后结合 GR，GR 进入细胞核，打开或关闭一批靶基因	FKBP5、TSC22D3 等就是经典的靶基因
细胞系 / 供者	不同人的细胞基线表达就不一样，这和处理无关	样本表的 `cell` 列
配对设计	同一个细胞系处理前后各测一次，比较的是"自己和自己"	公式 `~cell + dex`

二、统计概念

概念	为什么需要	在代码里
文库大小（library size）	样本 A 测了 3000 万条读段，样本 B 只测了 2000 万条，原始计数不能直接比	`counts.sum(axis=0)`
CPM	每百万读段里有多少条落在这个基因上，用来粗略地消除测序深度的差异	`counts / counts.sum() * 1e6`
log 变换	表达量跨了好几个数量级。取对数后，倍数变化就成了加减，画图、做 PCA 都更合适	`np.log2(cpm + 1)`
低表达过滤	计数只有 0～3 的基因，比例变化全是噪声，还会增加检验次数	`(counts >= 10).sum(axis=1) >= 4`
PCA	把两万维的样本压缩到两维来看，样本之间最大的差异在哪个方向	SVD
负二项分布	计数数据的方差比均值大（称为过度离散），DESeq2 用它给计数建模	PyDESeq2 内部
size factor	DESeq2 自己做的标准化，比 CPM 更稳健，不会被少数超高表达基因带偏	`dds.deseq2()` 内部
log2 fold change	处理组比对照组高了多少倍，取 log2。1 表示 2 倍，-1 表示一半	`log2FoldChange` 列
p 值	假如真的没有差异，看到这么大差异的概率	`pvalue` 列
多重检验 / FDR	做两万次检验，p < 0.05 的假阳性大约就有 1000 个。BH 方法控制的是"显著结果里假阳性所占的比例"	`padj` 列
火山图	横轴是效应大小（log2FC），纵轴是证据强度（-log10 padj），右上角和左上角就是要找的基因	`volcano.png`
MA 图	横轴是平均表达量，纵轴是 log2FC。用来检查低表达基因的倍数是不是被放大了	`ma.png`

三、Python 要点

知识点	例子	用在
读压缩 TSV	`pd.read_csv(p, sep="\t", compression="gzip")`	任务 1、2
逐元素映射	`series.map(func)`、`series.map(dict)`	任务 1
按索引对齐	`df[cols] / series`：Series 的索引会自动对齐到 df 的列	任务 2
字符串批量处理	`index.str.split(".").str[0]`、`index.str.endswith("_PAR_Y")`	任务 2
正则 findall	`re.findall(r'(\w+) "([^"]*)"', attrs)` 返回键值对列表	任务 2
布尔筛选	`counts[(counts >= 10).sum(axis=1) >= 4]`	任务 3
按列除法	`counts / counts.sum(axis=0)`（每列除以本列的和）	任务 3
SVD	`u, s, vt = np.linalg.svd(x, full_matrices=False)`	任务 4
组合条件	`(a < 0.05) & (b.abs() >= 1)`，要加括号	任务 5

参考实现：de

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

展开参考解

看答案

"""实验 2 · 批量 RNA-seq 差异表达（参考解）

数据：airway（Himes et al. 2014, PLoS One；GEO GSE52778），4 个人气道平滑肌细胞系，
每个细胞系一份未处理、一份地塞米松（dexamethasone，糖皮质激素）处理。
问题：地塞米松让哪些基因的表达发生了变化？
"""
from __future__ import annotations

import gzip
import re
import sys
import warnings
from pathlib import Path

import numpy as np
import pandas as pd

HERE = Path(__file__).resolve().parent
ROOT = HERE.parents[2]
sys.path.insert(0, str(ROOT / "tools"))
import datasets  # noqa: E402

# 只比较这两组；原研究还有沙丁胺醇（albuterol）处理组，本实验不用
TREATMENTS = {"Untreated": "untreated", "Dexamethasone": "dex"}

# 文献里反复验证过的糖皮质激素应答基因（地塞米松处理后上调），用来检验分析是否靠谱
KNOWN_DEX_UP = ["FKBP5", "TSC22D3", "PER1", "ZBTB16", "KLF15", "DUSP1"]


# ---------------------------------------------------------------- 任务 1：样本表
def parse_attributes(text: str) -> dict[str, str]:
    """recount3 的属性串 'cell line;;N61311|treatment;;Untreated' → {'cell line': 'N61311', ...}"""
    out = {}
    for item in text.split("|"):
        if not item:
            continue
        key, _, value = item.partition(";;")
        out[key.strip()] = value.strip()
    return out


def load_samples(sample_md_path, qc_md_path) -> pd.DataFrame:
    md = pd.read_csv(sample_md_path, sep="\t", compression="gzip")
    qc = pd.read_csv(qc_md_path, sep="\t", compression="gzip")
    attrs = md["sample_attributes"].map(parse_attributes)
    samples = pd.DataFrame({
        "run": md["external_id"],
        "cell": attrs.map(lambda a: a["cell line"]),
        "treatment": attrs.map(lambda a: a["treatment"]),
    })
    samples = samples[samples["treatment"].isin(TREATMENTS)].copy()
    samples["dex"] = samples["treatment"].map(TREATMENTS)
    avg_len = qc.set_index("external_id")["star.average_mapped_length"]
    samples["avg_len"] = samples["run"].map(avg_len)
    samples = samples.drop(columns="treatment").set_index("run")
    return samples.sort_values(["cell", "dex"], ascending=[True, False])


# ---------------------------------------------------------------- 任务 2：计数矩阵与注释
def coverage_to_counts(coverage: pd.DataFrame, avg_len: pd.Series) -> pd.DataFrame:
    """recount3 存的是"碱基覆盖度之和"，除以平均比对长度才约等于读数，再四舍五入成整数。"""
    return (coverage[avg_len.index] / avg_len).round().astype("int64")


_ATTR = re.compile(r'(\w+) "([^"]*)"')


def load_gene_annotation(gtf_path) -> pd.DataFrame:
    rows = []
    with gzip.open(gtf_path, "rt", encoding="utf-8") as f:
        for line in f:
            if line.startswith("#"):
                continue
            cols = line.rstrip("\n").split("\t")
            if cols[2] != "gene":
                continue
            attrs = dict(_ATTR.findall(cols[8]))
            rows.append((attrs["gene_id"], attrs.get("gene_name", ""), attrs.get("gene_type", ""),
                         cols[0], int(cols[5])))
    ann = pd.DataFrame(rows, columns=["gene_id_versioned", "gene_name", "gene_type", "chrom", "length"])
    ann = ann[~ann["gene_id_versioned"].str.endswith("_PAR_Y")]
    ann.index = ann["gene_id_versioned"].str.split(".").str[0].rename("gene_id")
    return ann.drop(columns="gene_id_versioned")


def load_airway() -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
    """读入并整理好：(counts 基因×样本, samples 样本表, annotation 基因注释)。"""
    samples = load_samples(datasets.path_of("airway_sample_md"), datasets.path_of("airway_qc_md"))
    coverage = pd.read_csv(datasets.path_of("airway_gene_sums"), sep="\t", comment="#",
                           index_col=0, compression="gzip")
    coverage = coverage[~coverage.index.str.endswith("_PAR_Y")]
    coverage.index = coverage.index.str.split(".").str[0].rename("gene_id")
    counts = coverage_to_counts(coverage, samples["avg_len"])
    ann = load_gene_annotation(datasets.path_of("gencode_genes"))
    return counts, samples, ann.reindex(counts.index)


# ---------------------------------------------------------------- 任务 3：过滤与标准化
def filter_low_counts(counts: pd.DataFrame, min_count: int = 10, min_samples: int = 4) -> pd.DataFrame:
    """至少 min_samples 个样本的计数 ≥ min_count 才保留（min_samples 取最小组的样本数）。"""
    keep = (counts >= min_count).sum(axis=1) >= min_samples
    return counts[keep]


def cpm(counts: pd.DataFrame) -> pd.DataFrame:
    return counts / counts.sum(axis=0) * 1e6


def log_cpm(counts: pd.DataFrame, prior: float = 1.0) -> pd.DataFrame:
    return np.log2(cpm(counts) + prior)


# ---------------------------------------------------------------- 任务 4：PCA
def pca(expr: pd.DataFrame, n_top: int = 500, n_components: int = 2) -> tuple[pd.DataFrame, np.ndarray]:
    """expr: 基因×样本（log 尺度）。取方差最大的 n_top 个基因，按基因中心化后做 SVD。
    返回 (样本×主成分 的得分表, 各主成分解释的方差比例)。"""
    top = expr.loc[expr.var(axis=1).sort_values(ascending=False).index[:n_top]]
    x = (top.T - top.T.mean(axis=0)).to_numpy()          # 样本 × 基因，每列均值为 0
    u, s, _vt = np.linalg.svd(x, full_matrices=False)
    scores = u[:, :n_components] * s[:n_components]
    ratio = (s ** 2 / np.sum(s ** 2))[:n_components]
    cols = [f"PC{i + 1}" for i in range(n_components)]
    return pd.DataFrame(scores, index=top.columns, columns=cols), ratio


# ---------------------------------------------------------------- 任务 5：DESeq2
def run_deseq(counts: pd.DataFrame, samples: pd.DataFrame) -> pd.DataFrame:
    """配对设计 ~ cell + dex：先扣掉细胞系之间的差异，再看地塞米松的效应。"""
    from pydeseq2.dds import DeseqDataSet
    from pydeseq2.default_inference import DefaultInference
    from pydeseq2.ds import DeseqStats

    inference = DefaultInference(n_cpus=1)
    dds = DeseqDataSet(
        counts=counts.T,
        metadata=samples[["cell", "dex"]],
        design="~cell + dex",
        inference=inference,
        quiet=True,
    )
    with warnings.catch_warnings():
        # 8 个样本、5 个参数，残差自由度只有 3，PyDESeq2 会提醒离散度先验估计不稳——已知情况
        warnings.filterwarnings("ignore", message="As the residual degrees of freedom")
        dds.deseq2()
    stats = DeseqStats(dds, contrast=["dex", "dex", "untreated"], inference=inference, quiet=True)
    stats.summary()
    return stats.results_df


def classify(results: pd.DataFrame, padj: float = 0.05, lfc: float = 1.0) -> pd.Series:
    sig = results["padj"].lt(padj) & results["log2FoldChange"].abs().ge(lfc)
    out = pd.Series("ns", index=results.index)
    out[sig & (results["log2FoldChange"] > 0)] = "up"
    out[sig & (results["log2FoldChange"] < 0)] = "down"
    return out


def naive_fold_change(counts: pd.DataFrame, samples: pd.DataFrame) -> pd.Series:
    """旧版实验的做法：处理组原始计数均值 ÷ 对照组原始计数均值。用来和 DESeq2 对比，看它错在哪。"""
    dex = counts[samples.index[samples["dex"] == "dex"]].mean(axis=1)
    ctl = counts[samples.index[samples["dex"] == "untreated"]].mean(axis=1)
    return dex / ctl


# ---------------------------------------------------------------- 任务 6：出图与报告
def main(out_dir) -> None:
    import matplotlib

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

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

    counts, samples, ann = load_airway()
    kept = filter_low_counts(counts)
    logexpr = log_cpm(kept)
    print(f"基因：{len(counts)} → 过滤低表达后 {len(kept)}；样本：{len(samples)}")

    # PCA
    scores, ratio = pca(logexpr)
    fig, ax = plt.subplots(figsize=(5.5, 4.5))
    sns.scatterplot(data=scores.join(samples), x="PC1", y="PC2", hue="dex", style="cell", s=90, ax=ax)
    ax.set(xlabel=f"PC1 ({ratio[0]:.0%})", ylabel=f"PC2 ({ratio[1]:.0%})", title="PCA of log2 CPM (top 500 var genes)")
    ax.legend(fontsize=7, loc="best")
    fig.tight_layout()
    fig.savefig(out / "pca.png", dpi=150)
    plt.close(fig)

    # DESeq2
    res = run_deseq(kept, samples)
    res = res.join(ann[["gene_name", "gene_type"]])
    res["class"] = classify(res)
    res.sort_values("padj").to_csv(out / "de_results.csv", encoding="utf-8-sig")
    n_up, n_down = (res["class"] == "up").sum(), (res["class"] == "down").sum()
    print(f"显著上调 {n_up}，下调 {n_down}（padj < 0.05 且 |log2FC| ≥ 1）")

    colors = {"up": "#d62728", "down": "#1f77b4", "ns": "#bbbbbb"}
    # 火山图
    fig, ax = plt.subplots(figsize=(6, 5))
    y = -np.log10(res["padj"].clip(lower=1e-300))
    ax.scatter(res["log2FoldChange"], y, c=res["class"].map(colors), s=4, linewidths=0)
    for g in KNOWN_DEX_UP:
        row = res[res["gene_name"] == g]
        if len(row):
            ax.annotate(g, (row["log2FoldChange"].iloc[0], y[row.index[0]]), fontsize=7)
    ax.axhline(-np.log10(0.05), ls="--", lw=0.6, color="gray")
    ax.axvline(1, ls="--", lw=0.6, color="gray")
    ax.axvline(-1, ls="--", lw=0.6, color="gray")
    ax.set(xlabel="log2 fold change (dex vs untreated)", ylabel="-log10 adjusted p", title="Volcano plot")
    fig.tight_layout()
    fig.savefig(out / "volcano.png", dpi=150)
    plt.close(fig)

    # MA 图
    fig, ax = plt.subplots(figsize=(6, 4.5))
    ax.scatter(np.log10(res["baseMean"] + 1), res["log2FoldChange"], c=res["class"].map(colors), s=4, linewidths=0)
    ax.axhline(0, color="black", lw=0.6)
    ax.set(xlabel="log10 mean normalized count", ylabel="log2 fold change", title="MA plot")
    fig.tight_layout()
    fig.savefig(out / "ma.png", dpi=150)
    plt.close(fig)

    # 前 30 个差异基因热图（每个基因做 z-score，看模式而不是看绝对值）
    top = res[res["class"] != "ns"].sort_values("padj").head(30)
    z = logexpr.loc[top.index]
    z = z.sub(z.mean(axis=1), axis=0).div(z.std(axis=1), axis=0)
    z.index = top["gene_name"].where(top["gene_name"] != "", top.index)
    z.columns = [f"{samples.loc[c, 'cell']}_{samples.loc[c, 'dex']}" for c in z.columns]
    grid = sns.clustermap(z, cmap="vlag", center=0, figsize=(7, 9), col_cluster=True, row_cluster=True)
    grid.savefig(out / "heatmap_top30.png", dpi=150)
    plt.close(grid.figure)

    # 旧方法 vs DESeq2
    naive = naive_fold_change(kept, samples)
    naive_hits = naive[naive > 2].index
    naive_ns = (res.loc[naive_hits, "padj"].fillna(1) >= 0.05).sum()

    known = res[res["gene_name"].isin(KNOWN_DEX_UP)][["gene_name", "baseMean", "log2FoldChange", "padj"]]
    lines = [
        "# 实验 2 结果摘要（airway：地塞米松 vs 未处理，配对设计 ~ cell + dex）\n",
        f"- 基因：{len(counts)} → 过滤后 {len(kept)}；样本 {len(samples)}（4 细胞系 × 2 处理）",
        f"- PCA：PC1 解释 {ratio[0]:.0%}，PC2 解释 {ratio[1]:.0%}",
        f"- 显著差异（padj < 0.05 且 |log2FC| ≥ 1）：上调 {n_up}，下调 {n_down}",
        f"- 旧方法（原始均值比 > 2 倍）挑出 {len(naive_hits)} 个基因，其中 {naive_ns} 个在 DESeq2 里并不显著\n",
        "## 已知糖皮质激素应答基因\n",
        "| gene | baseMean | log2FC | padj |",
        "|---|---:|---:|---:|",
        *[f"| {r.gene_name} | {r.baseMean:.0f} | {r.log2FoldChange:.2f} | {r.padj:.2e} |"
          for r in known.itertuples()],
    ]
    (out / "结果摘要.md").write_text("\n".join(lines) + "\n", encoding="utf-8")
    print(known.round(4).to_string(index=False))

← 实验 1：序列基础（FASTA · GC 含量 · 反向互补 · 翻译 · ORF）目录实验 3：富集分析（ORA 超几何检验 · BH 校正 · GSEA） →