从「未实测」到 217 个契约产物：一次不调参的生信复现与验收

这座档案室里的项目大多有两种来路：一种是从零开始建一个能用的东西，另一种是把别人说「可以这样」的事，真正做一遍看会发生什么。这一期属于后者——它的起点是两份文档，两份把自己的「没跑过」坦诚写在脸上的文档。

整理时间：2026-10-11 ｜ 项目：biolearn-practice（GitHub） 运行环境：Ubuntu 24.04.5 虚拟机 · 4 vCPU / 15.6 GiB 内存 / 48 GiB 磁盘 ｜ 全部隔离环境建在用户目录，系统 R/Python 零改动 作业窗口：2026-10-10 19:40 → 10-11 00:39 UTC（约五小时：环境安装 → 全部分析 → 验收 → 首次发布） 说明：本文所有数字取自项目实际产物、验收报告与 SHA 清单，非估算；参考文档数字仅供对照，不作为验收输入。

0. 起点：两份把「没跑过」写在脸上的文档

10 月 10 日，两份学习文档放到了桌上：一份是 5 周路线的简版，一份是逐日拆解的详细实践指南。它们规划了一个 RNA-seq + 单细胞的复现项目——用 airway 数据跑三种 bulk 差异分析工具，用 PBMC3k 数据跑两条单细胞路线，最后再用 nf-core 官方测试收尾。

专业的地方在于：这份指南给每一句话都挂了核实状态。

也就是说：图纸是完整的，但其中一半没有施工记录。文档还给出了参考数字——FDR<0.1 显著基因 4,381、QC 后 2,638 个细胞、Seurat 9 簇、Scanpy 8 簇——以及一句把它自己交给现实的作业：

「2,638 个细胞，8 个 Leiden 簇，写出 9 对 8 的原因分析」

于是这个项目要做的就不是「再写一份文档」，而是：在一台真实虚拟机上，把那行【未实测】逐条摘掉——并且不靠修改数字来完成。

1. 先立规矩：什么叫「做完」

动工之前先写约定，是这个项目的第一个动作。一组规则被写进 AGENTS.md，其中两条决定了后面所有故事的走向：

• 验收数字必须从本次实际产物读取，记录数据、参数、环境版本及产物路径；不得将参考数字写成实测结果。
• 不为匹配参考数字而调参；参数选择应有方法或数据依据，差异应如实记录。

其余条款同样明确：参考材料待核实、优先复用隔离环境不覆盖系统 R/Python、不自动推送 GitHub / 不改代理 / 不删旧资料、进度文档如实记录阻塞与下一步。

与约定配套的，是一套五档的证据分级，后面每一份报告都按它来写：

⚠️ 教训：这套分级看似只是措辞纪律，实际是整个项目的脊柱。 它提前堵死了所有「差不多算通过」的模糊地带——WARN 必须带着解释留着，NOT RUN 必须写明最小解决方案，不许含糊成 PASS。

2. 环境：三间实验室，都在用户目录里

复现的第一关是环境。原则一条：不碰系统，全部隔离。 Miniforge 装进用户目录，然后在它下面建起三间用途清晰的实验室，外加三个辅助环境：

每个环境都导出了带 build 与频道的 explicit 锁文件；每一轮实际运行时再导出一次「本次实际清单」，与锁文件逐项核对——本次验收中两边完全一致。系统 Python（3.12.3）原样保留，R 与 Conda 不进 PATH，一律用绝对路径调用。

环境审计同样是实测口径：核心包逐个 library() / import 验证（不用查文件式的假验证），七个网络目标逐个 HEAD 请求返回 200，TLS 校验通过；连 scikit-misc 的 conda 版本 0.5.3 与 Python 元数据 0.0.0 不一致这种小事，也单独记录为 WARN，而不是「手改一下元数据让它好看」。

环境不是复现的背景条件，它是实验记录的一部分。

3. 契约工程：把「运行」变成可验收的流程

分析能跑通只是第一步。这个项目真正的骨架，是把整个流程变成「先说好要交付什么，再让机器去满足它」。

Snakefile 只有四条规则，但配合两份契约文件，定义了全部战场：

• 主契约 workflow/outputs.json：bulk 69 + PBMC3k 97 = 166 个主产物
• 诊断契约 workflow/diagnostic-outputs.json：34 + 17 = 51 个补充产物
• 合计 217 个声明产物——少一个、空一个，契约检查都会拦下

输入（airway 包 / 10x PBMC3k / nf-core 固定测试数据）
   │
   ▼
Snakefile（4 条规则）
   │  契约：166 主产物 + 51 诊断产物 = 217 个声明产物
   ▼
每次运行 = 一个独立 attempt（.attempts/<路线>-<时间戳>/）
   │  逐阶段记录：命令 · 脚本 SHA-256 · 起止时间 · 退出码
   ▼
契约检查 + 内部验收 ──通过──▶ 发布（旧版本移入 history/，新版本进正式目录）
   │
   └──失败──▶ attempt 原样保留，不进入交付目录；已发布产物哈希不动

跑每一条路线的 run_matrix.py 是个很小的脚本，却把「可信」两个字拆成了动作：每次运行先建独立的 .attempts/ 沙盒；每个阶段跑完就把命令、脚本哈希、起止时间和退出码写进 run-provenance.json；任一阶段非零退出就立即终止，一个字节都不会发布；全部通过后才把结果发布进正式目录，旧版本整体移入 history/。

这套机制后来被两个实验正面检验过：

其一，失败隔离实测。 把已验证的 bulk 输出复制成测试 fixture，再把环境命令故意指向 /bin/false——工作流如期以 exit 1 失败，已发布的 75 个文件哈希全部未动，failed attempt 的 completed=false 如实落盘。这道证明的意思是：环境坏了，也不该污染证据。

其二，缺图重建证明。 把新生成的 bulk/pca.png 从产物目录移走（不是删除，原图保留在证据目录），重新运行工作流。Snakemake 精确判定只有 bulk 规则需要重跑，重建出的 PCA 图与移走的原图 SHA-256 完全一致（041bf3ae…b84ca9）；五张主要统计表逐字节相同，其余三条路线的 157 个文件哈希零变化。最终 dry-run 显示 Nothing to be done。

⚠️ 教训：「可复现」不是重跑一次看着差不多，而是删掉一个产物，机器能把它重建得与原来逐字节相同。 SHA 比「我验证过了」这四个字硬得多。

发布之后还有一层：验收器必须独立于分析器。 README 对照检查 16/16、扩大核验 158/158、补充诊断 18/18——而扩大核验器自己的第一版还翻过车：它把 raw RDS 的 list 误判成 Seurat 对象、把注释计划的参数也拿去当脚本比 SHA，跑出 156/158。修好验收器后重跑才是 158/158，两版结果与失败日志都保留在案。

⚠️ 教训：验收器本身也要被验收。 它的 bug 不是丑闻，是必须归档的证据——否则你无法区分「分析失败」和「量具失准」。

这一段收尾时的完整数字是：217 个声明产物全部存在、非空、可读；连同附加文件回读 232/232；三个 RDS、五个 H5AD 实际打开逐元素核对；历史 results 的 7,790 个文件 SHA 全部与基线一致——没有为了「干净」删掉任何旧东西。

4. bulk：同一个问题，三种回答

airway 数据来自人气道平滑肌细胞的 4 个供体、每个供体一对地塞米松处理/未处理样本，共 8 个样本、63,677 个基因。三个工具跑同一份计数矩阵、同一个设计 ~ cell + dex、同一条显著性线（FDR<0.1），但各自保留自己的方法学：

实际结果：

三个结果互相之间的对照，比单个数字更能说明问题：

谱系一目了然：DESeq2 与 PyDESeq2 是同一家族的两个实现，集合相似度 0.936；edgeR 走的是过滤与检验都不同的另一条路，重合度约 0.70–0.74——但三种方法给出的基因方向没有一个是反的。方法差异被如实保留，而不是被抹平成一句「结果一致」。

然后是这一期的核心剧情：参考数字是 4,381，实测是 4,916。差值 +535。

按学习路径文档的口径，4,381 来自 rnaseqGene 工作流（Salmon 量化 + tximeta），和本项目直接用 airway 现成计数矩阵不是同一条输入链路。项目对这笔差值的处理只有三条：重复跑一遍得到同样的 4,916；把差值原样写进 reference_diagnostics.json；以及一行 policy——

「diagnostic only; no pass/fail threshold and no parameter tuning」 （仅作诊断；不设通过阈值；不调参数）

没有阈值、没有调参、没有强行归因。附带的诊断里还有一个有意思的细节：参考文档列出的 6 个 top 基因（MAOA、DUSP1、SAMHD1、SPARCL1、SERPINA3、GPX3）中，实测 top20 命中了 5 个。数字对不上，但信号方向是对的——这正是「诊断」该有的用法：它提示你去看，而不是命令你去改。

另有一个 WARN 被完整保留：PyDESeq2 警告残差自由度「小于 3」。查源码后发现条件是 (num_samples - num_vars) <= 3，8−5=3 确实触发，属于小样本先验估计的正常局限——记录、解释、保留，但不通过删掉配对项来消警告。

5. 单细胞：2638 个细胞上的两条路线

PBMC3k 是 10x Genomics 的经典外周血单核细胞数据：2,700 个细胞、32,738 个基因。两条路线读同一份已过滤 UMI 矩阵，先做统一的严格 QC——检测基因数 200–2500 之间、线粒体比例低于 5%——最后都保留 2,638 个细胞、13,714 个基因（拒绝 62 个：5 个高基因数、57 个高线粒体）。按条形码精确对齐：共有 2,638 个，双方独有均为 0。

从完全相同的细胞出发，两条路线各自按教程跑完整流程：

结果：两条路线都得到 9 个簇——但分出来的不是同一张地图：

• ARI（调整兰德指数）= 0.6279，NMI = 0.7545
• 两边都是 9 簇，但交叉列联表显示主要分歧在 T 细胞区域（同一个细胞群被切成了不同边界）
• 两套推定的细胞类型注释一致率 0.9500——注意，这不是准确率，只是两套标签的吻合度；注释本身始终标注为 putative（推定）

高变基因集合也不一样：共有的只有 1,462 个。

⚠️ 教训：相同的簇数可以掩盖不同的结论。 「9 对 9」看起来比「9 对 8」整齐，但 0.628 的 ARI 说明这只是两个近似的划分——数字对齐从来不等于结论一致。这也是项目坚持不拿簇数当验收目标的原因：簇数是现象，不是目标。

两条路线的 marker 检验（cluster 内 vs 其余细胞的 Wilcoxon 检验）刻意不做直接对比：Seurat 返回 6,964 行、Scanpy 返回 123,426 行——检验全集与校正方式不同，行数差不能当作一致性指标。这条口径写在方法文档里，与「不为数字调参」同源。

6. 「9 对 8」：把一句作业变成一组实验

回到学习文档留下的那道作业：教程里 Seurat 得 9 簇、Scanpy 得 8 簇，请「写出 9 对 8 的原因分析」。

这道作业最常见的完成方式，是写一段「因为算法不同、参数不同」的解释。项目选了另一条路：把候选原因逐条做成对照实验，让结论向证据收敛。

先是一组控制变量实验——每次都只改一个旋钮：

然后是三组补充实验：Scanpy 邻居 k 从 10 调到 20 → 8 簇；把近似邻居换成精确计算（k=10）→ 仍然是 9 簇、ARI 1.0（这个因素被排除）；换用 Seurat 的高变基因集 → 8 簇。

最后是图诊断——把两条路线「底层地基」直接摊开比：Seurat 的 SNN 图与 Scanpy 的 fuzzy 图，公共边 Jaccard 只有 0.055（图本身差得很远）；但把同一个 Leiden 算法分别放到两张图上：分辨率 0.5 时 SNN 给 9 簇、fuzzy 给 6 簇；分辨率 1.0 时给 11 与 9 簇。其中 SNN + res.5 的划分与 Seurat 原结果的 ARI 高达 0.961。

这些实验给出的结论被谨慎地收在实测边界内：分辨率、种子、图公式、PC、回归、HVG 都可以影响簇数与划分；不能把全部差异归因于单一因素，也没有把 8 簇选来当「匹配参考」的新主结果。主参数保持原样。

⚠️ 教训：解释「为什么不同」，要用实验，不要用形容词。 而实验的结论要允许它是「多个因素交织、无法唯一归因」——诚实的边界比漂亮的答案更值钱。

顺带一提：实测结果是 9 对 9，参考文档写的 8 从未在本机出现过。项目也没有为了这句「9 对 8」把任何参数调成 8。

7. 计算能回答什么：doublet 与参考映射

主分析之外，项目做完了两份「诊断性」验证——它们回答的都不是「我对不对」，而是「一个方法能测到什么程度」。

doublet（双重细胞）检测。 Scrublet 在 2,700 个细胞上标记 36 个疑似双重细胞（主 QC 后剩 34 个），阈值 0.2460。项目对这份结果的表述极其克制：检测比例不是实验 doublet 率（期望率 0.05 是假设值而非真值）；没有任何细胞被据此删除，只作为诊断记录到每个簇的分布表里。

独立参考映射。 用 Scanpy 随包的 PBMC68k reduced 参考（700 细胞 / 765 基因，与查询共有 764 基因）做距离加权 kNN 映射，留出 30%（210 细胞）做验证：holdout accuracy = 0.8429，macro-F1 = 0.8348；两条主路线与参考的标签一致率分别为 0.813（Seurat）与 0.805（Scanpy）。

但这份验证的 WARN 同样被醒目地保留：该参考缺失 CD14 / MS4A7 / PPBP / PF4 这几类关键 marker，没有 FCGR3A 特异的单核细胞类——所以它只能提示、不能判卷。项目最后一行写着：「计算参考不能替代实验真值」，蛋白验证 NOT RUN（缺数据）。

⚠️ 教训：知道一个方法能回答什么，和知道它回答不了什么，同等重要。 把「计算上与参考 80% 一致」说成「注释被验证了」，是这个领域最常见的越界。

8. nf-core：一次被归档的失败，和 234/234 的通过

矩阵分析之外，文档还要求跑一次 nf-core/rnaseq 官方测试——真实的 Nextflow 流程、真实的容器、真实的多工具编排。版本被钉死：rnaseq 3.27.0，commit a1fcdddd，test+docker 配置。

第一次尝试，失败，且失败得很标准： 流程退出码 1，首个任务 GUNZIP_GTF 退出 127，报错原文是 docker: command not found。整份记录被原样归档：completion-status 标记 blocked_container_runtime_missing、MultiQC 未生成，目录命名为 attempt-1 保留。没有假装它没发生过。

容器引擎需要有人到虚拟机终端里装——这一步由人在终端里完成，凭据从未经过对话。Docker 装好后新的问题立刻出现：daemon 直接拉镜像被拒绝（dial tcp …: connect: connection refused）。项目没有去改系统代理，而是请出一个独立的 Skopeo（1.24.1）加上项目自带客户端，继承已有网络环境完成精确镜像导入，再用一个 Alpine 容器做预检：容器可运行、DNS 正常、容器内 HTTPS 请求官方样本表返回 200。

第二次尝试，全绿：

收尾时还有一关：MultiQC 官方报告里的图表颜色异常。根因定位到展示层——MultiQC 1.35 把裸 RGB 字符串（如 55,126,184）交给图表库，93 次转换警告后回退成黑色。修复方式被严格限定为「只修展示」：写一个独立 adapter 把裸 RGB 转成等值的 rgb(...) 格式，用原始任务的输入重新生成一份展示报告；36 个数据模块与 14 组统计与官方 JSON 完全相等，科学数据一个字节没动。 前三次失败的展示尝试（目录缺失、容器环境变量、旧报告被重复扫描）也全部保留，只用第 4 次有效结果。

⚠️ 教训：展示层的问题只修展示层。 数据与呈现的分界线，一旦模糊，「修复」就会变成「篡改」。此外还要记住一句底线——官方小型测试 ≠ airway 原始 FASTQ 的完整复现，这句话在每个相关文档里都被写死在结论旁。

9. 发布：Git 收什么、不收什么

最后一小时，项目上了 GitHub。发布边界本身就是一个决策现场：提交代码、文档、环境锁、关键图表与验收记录；不提交下载数据、隔离环境本体、RDS/H5AD、BAM/FASTQ、流程缓存、失败尝试与事务沙盒。提交前对候选集做了实测审计：约 64.8 MB、1,488 个文件，最大单文件约 13.2 MB，未发现密钥或典型令牌。

仓库公开在 github.com/Andrew-li777/biolearn-practice。与所有复现项目一样，一句边界被明确写在记录里：Git 发布不等于跨主机新环境复现——被排除的大对象（原始数据、序列化对象、完整历史结果）仍需另行归档，而在第二台机器上恢复锁环境重跑，至今仍是 UNVERIFIED。

10. 复盘：十条被验证的原则

这一夜的所有关键决策，收敛成十条：

1. 参考数字是诊断，不是分数线。 4,381 与 4,916 的差值被记录、被解释边界，但从不驱动参数——policy 原文只有一行：「仅作诊断；不设通过阈值；不调参数」。

2. 「不调参」要写成纪律，不能靠自觉。 它被写进项目约定，与「验收数字必须从实际产物读取」并列。凡是没写下来的自律，在数字诱惑面前都会失效。

3. 契约先行。 先声明 217 个产物，再让流程满足它；产物清单就是验收边界，缺一个、空一个都被拦下。「做完」因此有了一个不会含糊的定义。

4. 失败要归档，但不要留在交付目录。 attempt 沙盒里保留着每一次翻车：156/158 的验收器 bug、nf-core 首次退出 1、三次展示修复失败……交付目录里只有通过验收的版本。失败是证据，不是垃圾；但它也确实不是交付物。

5. 验收器本身也要被验收。 独立的扩大核验跑了 158 项检查——而它自己的第一版就是被修过两处误判才上岗的。量具失准与分析失败，是两种必须分开的事故。

6. SHA 比日志更硬。 删掉 pca.png，工作流重建出与原来逐字节相同的图；7,790 个历史文件以哈希证明未被触碰。「可复现」的终极形态是数学，不是回忆。

7. 相同数字可以掩盖不同结论。 两条单细胞路线都是 2,638 细胞、9 个簇，ARI 却只有 0.628。对齐的数字会撒谎，对齐的方法不会。

8. 计算参考能提示、不能定义真值。 参考映射约 80% 一致、缺了关键细胞类，它的正确定位是「诊断」，而不是「判卷人」；实验真值没有到位之前，注释永远标着 putative。

9. 展示层的问题只修展示层。 MultiQC 的颜色缺陷用独立 adapter 解决，数据零改动、原始报告零覆盖——科学结果与呈现方式之间，需要一条防火墙。

10. 边界写出来才存在。 PASS 之外，WARN / NOT RUN / UNVERIFIED 各有定义与清单。一个项目对「自己不知道什么」的表述精度，决定了它有多可信。

11. 现状与下一步

今天的样子：一个可在本机一条命令复现的矩阵流水线（snakemake --cores 4 --resources mem_mb=14000 --config output_root=…）+ 一条包括官方流程测试在内的完整证据链。217 个声明产物、16/16 对照、158/158 扩大核验、18/18 补充诊断、267/267 官方测试核验；GitHub 仓库已公开发布。历史资料零删除，系统环境零改动。

按依赖排出的下一步：

1. 跨主机新环境复现（当前 UNVERIFIED）——在另一台授权主机上恢复锁环境、执行同一条命令、比较产物。锁文件、额外依赖与命令都已保存，缺的只是第二台机器。
2. 实验真值——注释的准确性需要蛋白或其他实验测量数据来对齐；在拿到数据之前，任何计算验证都替代不了它。
3. 容量与长期维护——大运行前复核磁盘水位；输入或脚本变化时由 Snakemake 依赖自动触发重跑。

复现的终点不是得到与参考一致的数字，而是得到一组能经受复核的证据。数字可以不同——证据链必须完整。

参考链接

• 项目仓库：https://github.com/Andrew-li777/biolearn-practice（代码 / 文档 / 环境锁 / 验收结果）
• 数据来源：Bioconductor airway 1.26.0（GSE52778）；10x Genomics PBMC3k 已过滤矩阵；nf-core/rnaseq 3.27.0 官方测试
• 方法文档（仓库内）：docs/bulk-methods.md、docs/pbmc3k-methods.md、docs/verification-report.md、docs/nfcore-test.md
• 教程依据：Seurat PBMC3k 教程、Scanpy clustering-2017 教程、Bioconductor rnaseqGene 工作流
• 姊妹档案：本刊 第 004 期 · 翻译文库——它完成的第一个 200 块整书项目，正是一本生信教材

本文记录 biolearn-practice 的复现与验收过程，所有数字取自实际产物、验收报告与 SHA 清单；参考文档数字仅作对照。跨主机复现与实验蛋白验证等项标记为未验证 / 未执行；完整清单以仓库内 docs/verification-report.md 的当前版本为准。