同一疾病多个独立公开 RNA-seq 队列的差异表达可重现性：以「同队列子抽样 vs 跨队列复制」的差距量化小样本转录组结论的乐观偏倚

1 · 研究问题

对同一疾病，用统一管线在 5 个及以上相互独立的公开 bulk RNA-seq 队列上做差异表达（differential expression, DE）分析，跨队列的显著基因列表重叠率，比在单一大队列内部做同等样本量子抽样（subsampling）所得的队列内重叠率低多少？这个差距是否随样本量 N 单调收敛，能否给出一条"要使跨队列重叠率达到 0.5 所需的最小 N"的经验曲线？

2 · 研究背景与空白

背景。 bulk RNA-seq 的差异表达分析是生物医学论文最常见的产出形式之一：比较病例组与对照组，报出一份"显著上调/下调基因列表"，再做富集分析（enrichment analysis）。这份列表的可信度取决于统计检验力（statistical power），而检验力由样本量、效应量（effect size）和基因表达的生物学变异共同决定。人类组织样本的个体间变异远大于细胞系实验，因此小队列（每组 3–10 例）的 DE 列表是否稳定，直接决定了大量已发表结论是否站得住。

已有工作到哪一步。 Gerstner 等人（PLOS Computational Biology, 2025，PMC12077797）用 18 个真实数据集做了约 18,000 次子抽样实验：从一个大数据集中反复抽出规模为 N 的小队列，两两比较其 DE 结果的一致性。结论是每组少于 5 例时可重现性极低，但"低可重现性不等于低精确率"——18 个数据集中有 10 个在 N>5 时中位精确率仍然较高。他们还给出了一套自助法（bootstrap）程序，用来从手头数据估计预期可重现性。另一条线是 Li 等人（Genome Biology, 2022）指出 DESeq2/edgeR 在人群样本上会产生远高于名义水平的假阳性。

空白在于：现有的可重现性估计几乎全部来自同一个数据集内部的子抽样。 同一队列的样本共享招募标准、组织取材流程、文库制备批次和测序平台，因此队列内子抽样系统性地低估了真实世界"另一个实验室重做一遍"的不一致程度。真正决定文献可信度的是跨队列复制率，而这个量至今没有以"同一疾病、统一管线、多队列"的形式被系统标定过。这个空白适合一年期学生课题：数据完全公开且已统一处理（无需比对原始 fastq）、方法学完全公开、结论正负都成立——若跨队列率接近队列内率，说明现有小样本文献比想象中稳健；若显著更低，则给出一个可直接引用的乐观偏倚系数。

3 · 可检验假设

• H1：在相同的每组样本量 N（N ∈ {3, 5, 8, 12, 20}）下，跨队列 DE 列表的 Jaccard 重叠率显著低于同队列子抽样重叠率，差值在 N=5 时 ≥ 0.10（绝对值），且 95% 自助法置信区间不跨 0。
• H2（机制假设）：这一差距的主要来源是队列层面的技术与人群构成异质性而非随机噪声——若对每个队列先做 limma::removeBatchEffect 之外的队列内标准化 + 秩变换，再在秩空间比较，差距应缩小 ≥ 40%；若缩小不足 20%，则说明差异来自真实的人群生物学异质性，而非可校正的批次效应（batch effect）。

4 · 量化验收标准

1. 方法学校验（硬门槛）：选定 3 篇公开了完整差异基因表的原始论文（其数据在本项目队列集合内），用自建管线在原文全样本、原文所述统计设定下重跑，要求本项目得到的显著基因列表与原文发表列表的重叠率 ≥ 70%（以原文列表为分母，FDR<0.05、|log2FC|>1 的口径对齐），且 log2 倍数变化（log2 fold change）的 Spearman 相关 ≥ 0.85。任一队列不达标即视为管线未通过，该队列剔除；若 3 篇中 ≥ 2 篇不达标，全部后续结论无效，转入 go/no-go 的降级路径。
2. 统计口径预先写死：所有 DE 检验一律报 Benjamini–Hochberg FDR，不报裸 p 值；主分析阈值固定为 FDR < 0.05 且 |log2FC| > 1，并另做 FDR < 0.1 与仅 FDR 无 FC 阈值两组敏感性分析。重叠率同时报 Jaccard 指数与"前 200 基因的秩相关"两个口径（前者对列表长度敏感，后者不敏感），两者都必须报。
3. 批次效应必须显式检查：每个队列先做主成分分析（PCA）与 sva::num.sv 估计潜在因子数，报告病例/对照分组在前 3 个主成分上的解释比例；若某队列的病例/对照与测序批次完全混杂（confounded，即批次可完全预测分组），该队列标注为"不可去混杂"并单独列出，不进入主分析。
4. 检验力与效应量预先估计：用 RNASeqPower 或基于队列内离散度（dispersion）的模拟，给出每个 N 下检测 log2FC=1、基线表达中位数的基因所需的检验力，写入结果表。样本量不足以支撑的结论一律标注。
5. 规模下限：至少 5 个独立队列，每个队列病例与对照各 ≥ 15 例（保证可抽到 N=12）；每个 (队列, N) 组合做 ≥ 100 次子抽样重复；总计 ≥ 3,000 次 DE 运行。
6. 可复现性：全部脚本（R + Snakemake 或纯 R 驱动脚本）、随机种子、会话信息（sessionInfo()）与中间结果的校验和开源于 GitHub，第三方 git clone 后可一键重跑降规模版本（≤ 4 小时）。

5 · 数据与工具

用途	来源 / 工具
统一处理的基因水平计数矩阵（首选，避开比对算力墙）	recount3（https://rna.recount.bio/ ，R 包 recount3）。覆盖 SRA/GTEx/TCGA 共约 75 万条人鼠 RNA-seq run，全部用 Monorail 统一处理到 Gencode v26 基因水平。单个研究的 RangedSummarizedExperiment 约 60k 基因 × 样本数，300 样本时对象约 150–250 MB，16 GB 内存充裕
备用统一计数源	ARCHS4（https://maayanlab.cloud/archs4/ ，HDF5 单文件，人类版约 30–40 GB，需核实当前版本体积；可用 rhdf5 按样本切片读取，不必整文件载入）
原始 GEO 队列与元数据	GEO（https://www.ncbi.nlm.nih.gov/geo/ ），用 GEOquery 抓取 GSE 的样本表型表。候选疾病与队列须在第 3 周锁定；候选之一为特发性肺纤维化（IPF）肺组织，已知候选 GSE150910、GSE134692、GSE52463——三者的确切样本量与是否提供计数矩阵需核实
受控访问核查	recount3 中 TCGA 与 GTEx 的基因水平计数为开放数据；dbGaP 中的个体水平基因型、BAM/CRAM 与部分 GTEx 层级为受控访问，中学生无法申请，本项目一律不使用。所有选入队列必须能匿名下载，无需 DAC 审批
差异表达	DESeq2、edgeR（quasi-likelihood）、limma-voom。主分析用 limma-voom（单次运行秒级），DESeq2 用于校验与部分复核
批次与潜变量	sva、RUVSeq、limma::removeBatchEffect
富集分析（次要终点）	clusterProfiler + MSigDB Hallmark（50 条基因集，规模可控）
对照基准	Gerstner 等 2025 公布的队列内可重现性曲线；各原始论文发表的差异基因表。仅用于校验与对比，不计入本项目的数据贡献
算力口径	纯 CPU。limma-voom 单次（60k 基因 × 24 样本）< 5 s；DESeq2 单次 20–60 s。3,000 次 limma-voom ≈ 4 小时；若全用 DESeq2 ≈ 25–50 小时（可后台分夜跑）。从 fastq 重头比对为超出范围：300 样本 × 约 5 GB fastq ≈ 1.5 TB 下载，STAR 人类索引需约 32 GB 内存，HISAT2 单样本 4 核约 30–60 分钟，合计 150–300 机时——16 GB 笔记本不可行，方案中不设此路径

6 · 方法路径

1. 装环境（R 4.4 + Bioconductor），跑通 recount3 与 DESeq2 内置示例，先完成验收标准第 1 条的三篇论文复现，把重叠率与相关系数写入 validation/ 目录。
2. 按预先写死的筛选规则锁定疾病与 5–7 个队列（规则：同一组织、同为 polyA 或同为 rRNA-depleted 文库、病例/对照各 ≥ 15、计数矩阵可从 recount3 或 GEO 补充文件直接取得）。规则写在方案里，不得看到结果后再改。
3. 核实每个队列的批次-分组混杂情况与平台差异，出具 PCA 图与 num.sv 表；剔除不可去混杂队列，剔除理由逐条记录。
4. 建立两条并行的重抽样实验：(a) 队列内子抽样——在最大的单个队列内抽 N 例/组，重复 100 次；(b) 跨队列复制——在队列 i 抽 N 例/组、在队列 j 抽 N 例/组，对所有队列对重复 100 次。两条用完全相同的 DE 设定。
5. 计算每对结果的 Jaccard 指数与前 200 基因秩相关，用自助法给 95% 置信区间，画出重叠率–N 曲线（两条并列），并拟合达到 0.5 所需的 N。
6. 检验 H2：把队列内秩变换/批次校正加入管线，重跑跨队列分支，量化差距缩小幅度。
7. 独立交叉校验：换一个完全不同的统计框架（edgeR quasi-likelihood 与非参数的 Wilcoxon 秩和 + BH）重算主曲线，验证结论不依赖于 limma-voom 的模型假设。

7 · 新颖性边界

本课题不声称提出新的差异表达方法，不声称发现新的疾病生物学，不把任何单个队列的差异基因列表当作本项目的发现。"小样本 RNA-seq 可重现性低"这一定性结论已由 Gerstner 等（2025）、Li 等（2022）发表，学生不得声称为自己的发现。

已有工作具体完成了什么：Gerstner 等在 18 个数据集上做了约 18,000 次同一数据集内部的子抽样配对比较，得到队列内可重现性随 N 的曲线与一套自助法估计程序；他们的比较对象始终是同源样本。

本项目的贡献（且是主结论）：把评价维度从"队列内可重现性"换成"跨独立队列复制率"，并首次给出同一疾病下这两条曲线的并列差距，即量化"同队列子抽样对真实复制率的高估幅度"。这与研究手册里"随机划分 vs 时间划分对照量化高估幅度"是同一类贡献。

为什么有价值：读者引用小样本 RNA-seq 结论时，真正需要的是"另一个实验室能否重复"，而现有可重现性估计给的是同源样本的上界。给出这个差距系数，等于给现有文献配一把折扣尺。

风险声明："两条曲线差异不显著"同样是有效结论，但必须给出自助法误差棒，证明本设计在 N=5 时有能力分辨 0.10 的重叠率差异（预实验需先做这个功效检查）。

8 · 决策门槛（go / no-go）

• 第 4 周末：确认能否找到满足筛选规则的 ≥ 5 个独立队列。若合格队列 < 5，立即降级：其一，把疾病换成公开队列最多的方向（结直肠癌肿瘤/癌旁配对，recount3 中 SRA 项目数量级更大）；其二，把"独立队列"的定义放宽为"同一疾病的不同组织采样部位"，并在论文中显式标注这一降级及其对结论解释的限制。两条降级路径都保留"两条曲线并列 + 差距量化"的主结论框架。
• 第 8 周末（硬门槛）：验收标准第 1 条必须通过。若 3 篇复现中 ≥ 2 篇重叠率 < 70%，立即降级：改用作者已公开完整分析代码的队列（如带 GitHub 仓库的论文），把校验目标从"复现论文结果"改为"复现论文代码的输出"，阈值提高到 ≥ 90%。主结论框架不变。
• 第 20 周末：若跨队列曲线的自助法置信区间宽度 > 0.20（即无法分辨假设中的 0.10 差异），立即降级：把队列对数从"所有两两组合"扩到"含重复抽样的更大配对池"，或把 N 网格收缩到 {5, 12, 20} 三点以把重复次数从 100 提到 300。
• 第 36 周：结果冻结，进入英文写作与查重。第 44 周完成英文 PPT 初稿。
• 需提前核实而非边做边发现：(a) 候选 GSE 是否提供基因计数矩阵而非仅提供受控访问的原始数据；(b) recount3 对目标 SRA 项目的覆盖情况（并非所有 GEO 研究都在 recount3 中）；(c) 各队列的文库类型是否可比。三项必须在第 3 周前查清。
• 选择前提：适合已能独立写 R 脚本、且愿意接受"结论可能是'差距不显著'"的学生。本路线几乎不可能完全失败——即使 H1 与 H2 都被否证，两条并列曲线本身就是可发表的交付物。
• 预算裁剪顺序：先砍富集分析分支（次要终点），再砍 edgeR 交叉校验（保留 Wilcoxon），最后砍 N 网格点数。砍到只剩 limma-voom + 3 个 N 点时，主结论仍成立。

---