B03
单基因条形码物种界定方法间的不一致度能否被数据集属性预测：跨 40 个公开 COI 数据集的元分析与一条"何时可信"的可操作判据
1 · 研究问题
在 40 个及以上从 BOLD/GenBank 公开获取的 COI（细胞色素 c 氧化酶亚基 I）条形码数据集上，用统一管线跑四种主流物种界定（species delimitation）方法所得的分子操作分类单元（MOTU）数量之间的不一致度，能否由数据集本身在分析前就可测得的属性（每形态种序列数、序列长度覆盖、地理跨度、种内/种间距离间隙的宽度）线性预测？能否据此给出一条"当属性 X 低于阈值时四方法必然分歧、不应据单基因下分类学结论"的可操作判据？
2 · 研究背景与空白
背景。 DNA 条形码用一段标准化基因（动物用 COI 约 658 bp）把未知个体归入物种。当形态学信息不足时，研究者用物种界定算法直接从序列推断"有几个物种"。主流方法分两类：基于遗传距离的 ABGD / ASAP（Assemble Species by Automatic Partitioning）、BIN，以及基于系统树的 GMYC（Generalized Mixed Yule Coalescent）、PTP/bPTP/mPTP。它们的理论前提不同，因此在同一数据上给出不同答案。
已有工作到哪一步。 大量单类群论文报告了这种分歧的量级：中国螳螂目（2025）用 ASAP / jMOTU / bPTP / GMYC 分别得到 32 / 58 / 68 / 60 个 MOTU；螽斯属 Gampsocleis 分别得到 6 / 13 / 10 / 23 个；斯里兰卡 Sericini 金龟（2022）明确报告"多种界定方法彼此拟合很差，且与形态种拟合也差"；而石首鱼科 Stelliferinae（2023）中 ABGD/GMYC/bPTP 却高度一致（30–32 个 MOTU，多数与有效种吻合）。也就是说，分歧有时极大、有时极小，这一点在文献里是公开的。
空白在于：没有人把这几十个各自独立的单类群结果放到统一管线下重做，并检验"分歧大小是否可由数据集属性预先预测"。 现有文献的结论停在"结果依类群而异，需结合形态"，这是一句无法执行的建议。研究者真正需要的是：拿到一个新数据集时，能否在跑方法之前就判断结果可不可信。这个空白适合学生课题：数据全部匿名公开、每个数据集算力在分钟级、方法完全公开、结论正负都成立（"不可预测"本身就是对该领域的重要警示）。
3 · 可检验假设
- H1：以"四方法 MOTU 数的变异系数（coefficient of variation, CV）"为不一致度指标，用数据集属性做多元回归，跨 40 个数据集的留一交叉验证 R² ≥ 0.35；其中"每形态种平均序列数"与"种内/种间距离间隙宽度"两个变量的回归系数在 BH 校正后 FDR < 0.05。
- H2（机制/备择假设）：不一致主要来自采样密度而非真实的分类学复杂度——若对每个数据集把每形态种的序列数随机稀释到 3 条，四方法的 CV 应普遍上升且跨数据集差异缩小（组间方差下降 ≥ 30%）。若稀释后组间差异不缩小，说明不一致由类群固有的进化速率异质性驱动，采样量再大也无法消除，这是更强的结论。
4 · 量化验收标准
- 方法学校验（硬门槛）：从上述文献中选 3 个已公开完整序列登录号与 MOTU 结果的数据集（候选：中国螳螂目、Gampsocleis、Stelliferinae——三者补充材料是否给出可复原的登录号列表需核实），用自建管线在原文所述参数下重跑，要求每个方法的 MOTU 数与原文报告值的相对偏差 ≤ 15%，且系统发育树上原文强调的关键分支（原文标注支持率 ≥ 0.95 的节点）在本项目重建树中的超快自展（UFBoot）支持率 ≥ 90%。3 个数据集中 ≥ 2 个不达标，则管线不可信，后续全部元分析结论无效。
- 统计口径预先写死：回归的显著性一律报 BH 校正后的 FDR，不报裸 p 值（40 个数据集 × 多个候选属性构成多重检验）。R² 必须报留一交叉验证 R²而非拟合 R²，并明确写出两者差值（过拟合幅度）。所有数据集级统计量给 1,000 次自助法（按数据集整体重抽样）95% 置信区间。样本量仅 40，因此在开工前先做效应量与检验力估计：在 n=40、5 个预测变量下，能以 80% 检验力检出的最小 R² 是多少（预期约 0.30，须自行用模拟核算并写入方案）；若最小可检出 R² 高于 H1 阈值，必须先把数据集数扩到 60。
- 分类学标签的地位：BOLD/GenBank 中的形态种名仅作为对照标签使用，不计入本项目的数据贡献；且必须报告标签本身的不确定性（BOLD 记录的鉴定者与鉴定等级字段），把"无鉴定者信息"的记录比例写入每个数据集的元数据表。
- 规模下限：≥ 40 个数据集，每个数据集 ≥ 8 个形态种且 ≥ 60 条序列；覆盖 ≥ 4 个动物门/纲，避免结论只适用于昆虫。
- 主观判断的可复现化：数据集的纳入/排除规则、序列修剪规则（比对两端缺失比例阈值）、异常序列（可能的 NUMT 或污染）的剔除判据，必须全部写成可执行的数值规则并随代码发布，不得人工逐条挑选。
- 可复现性：Snakemake 或 Makefile 驱动的完整管线、全部登录号清单、随机种子、各工具版本号开源；提供一个 5 数据集的降规模版本，第三方在 ≤ 2 小时内可重跑。
5 · 数据与工具
| 用途 | 来源 / 工具 |
|---|---|
| COI 条形码序列与形态种标签 | BOLD Systems 公开数据门户（https://www.boldsystems.org/ ），支持按分类群批量下载 TSV/FASTA（含 BIN、鉴定者、采集经纬度）。BOLD v5 的批量下载接口与配额需核实；备用入口为 GenBank/NCBI Nucleotide 按 COI[Gene] AND <taxon>[Organism] 检索 + efetch |
| 数据集属性中的地理跨度 | BOLD 记录自带经纬度字段；无坐标记录单独统计比例 |
| 序列比对 | MAFFT（--auto）。1,000 条 × 658 bp 约 10–30 s；3,000 条约 2–5 min。COI 为编码基因，另用 MACSE 或按密码子比对做一次一致性检查 |
| 遗传距离 | R 包 ape（K2P 与 p-distance 两种，均报） |
| 距离法界定 | ASAP（网页服务与本地二进制均有，1,000 条秒级）、ABGD（同级） |
| 树法界定 | IQ-TREE 2 + ModelFinder + 1,000 次 UFBoot。500 条 × 658 bp、4 线程约 2–8 min；3,000 条约 1–3 小时。mPTP（比 bPTP 快 1–2 个数量级，1,000 条分钟级）；GMYC 用 R 包 splits，需超度量树（ultrametric tree） |
| 超度量化 | ape::chronos（惩罚似然，秒–分钟级）。BEAST2 严格钟属超出范围：300 条序列、10⁷ 步 MCMC 在 4 核 CPU 上通常需 1–3 天/数据集，40 个数据集不可行；方案中改用 chronos 并在论文中明确标注这一近似及其对 GMYC 结果的可能影响（须用 3 个小数据集做 chronos vs BEAST2 的对照，量化 MOTU 数差异） |
| 元分析统计 | R：lm / glmnet（含正则化以防 n=40 过拟合）、boot、p.adjust(method="BH") |
| 对照基准 | 上述四篇已发表论文的 MOTU 计数与树拓扑。仅用于校验与对比，不计入本项目的数据贡献 |
| 伦理 | 全部为已公开的序列数据库记录，无人体数据、无活体动物操作、无需伦理审批 |
6 · 方法路径
- 装环境（MAFFT / IQ-TREE 2 / mPTP / ASAP / R+ape+splits），跑通各工具自带示例，先完成验收标准第 1 条的三数据集复现，把 MOTU 数与关键分支支持率的复现表写入
validation/。 - 按预先写死的数值规则从 BOLD 抓取候选数据集（≥ 8 形态种、≥ 60 序列、≥ 500 bp 有效长度），生成 ≥ 40 个数据集的清单与元数据表；纳入/排除全过程留日志。
- 对每个数据集计算 6–8 个分析前可得的属性：每形态种平均与最小序列数、有效比对长度、缺失率、地理跨度（最大成对大圆距离）、种内最大 K2P 距离的中位数、种内-种间距离间隙宽度（barcode gap）、形态种数、单序列种（singleton）比例。
- 统一管线跑四种界定方法，记录每数据集的四个 MOTU 数、CV，以及与形态种划分的匹配指标（调整兰德指数 ARI 与匹配比例）。
- 核实工具能力边界：在 3 个小数据集上做 chronos vs BEAST2 超度量树的 GMYC 结果对照（BEAST2 只跑这 3 个，成本可控），量化近似带来的 MOTU 数偏差并写入限制条款。
- 做元回归（属性 → CV），报留一交叉验证 R²、系数与 BH-FDR、自助法区间；按 H1 阈值判定，并把回归转成一条阈值判据（例如"每形态种序列数 < 4 且 barcode gap < 2% 时，四方法 CV 中位数 > 0.4，不应据单基因下结论"）。
- 检验 H2 的稀释实验；最后做独立交叉校验——用一套完全不参与建模的 10 个新数据集做外部验证，报判据在外部集上的命中率。
7 · 新颖性边界
本课题不提出新的物种界定算法，不声称发现任何新物种，不对任何类群做分类学修订，不声称"不同界定方法结果不一致"为本项目发现——这一点已由斯里兰卡 Sericini（2022）、中国螳螂目（2025）等多篇论文明确报告。丘奖 2021 年金奖论文《A molecular phylogeny of cavernicolous Oniscidea…》做的是具体类群的系统发育重建与新种描述，本课题与之路径完全不同：本项目不描述新种、不做形态学工作，交付物是跨类群的方法学判据。
已有工作具体完成了什么：各单类群论文在自己的数据上并列跑 2–5 种方法，报出各自的 MOTU 数与和形态种的吻合度，结论止于"依类群而异，建议结合形态学"。它们没有统一管线、没有跨数据集比较、没有把不一致度当成因变量去建模。
本项目的贡献（且是主结论）：把评价维度从"某类群应界定出几个物种"换成"方法间不一致度本身是否可预测"，并交付一条可在分析前使用的定量判据 + 外部验证命中率。这属于研究手册所列的"评价维度转换 + 方法层面的可复现性贡献"，在计算生物学中是被承认的贡献类型，但定位必须讲清楚，否则会被误读为重复劳动。
为什么有价值：条形码数据每年以百万条增长，绝大多数新数据集不会有配套形态学工作。一条"什么时候单基因结论不可信"的先验判据，直接影响这些数据的使用方式。
风险声明："属性无法预测不一致度"（R² 接近 0）同样是有效结论，它意味着不一致源于不可事前观测的因素，从而否定"用简单指标筛数据集"的做法。但必须由第 2 条的检验力估计证明本设计（n=40）有能力检出 R²=0.35，否则该结论无效。
8 · 决策门槛（go / no-go）
- 第 4 周末：确认能否从 BOLD/GenBank 稳定批量抓到 ≥ 40 个合格数据集。若合格数据集 < 25，立即降级：其一，放宽形态种下限到 6 并把类群范围扩到植物 rbcL/matK 与真菌 ITS（代价是需为每个标记重设距离阈值，须在论文中分标记报结果）；其二，把因变量从"四方法 CV"改为"两方法（ASAP vs mPTP）的 ARI 不一致度"，此时可用数据集数量大幅增加。两条降级都保留"属性 → 不一致度回归 + 外部验证判据"的主结论框架。
- 第 8 周末（硬门槛）：验收标准第 1 条必须通过。若 ≥ 2 个复现数据集的 MOTU 数偏差 > 15%，先排查比对与修剪参数；四周内仍不通过则降级：把校验目标改为"复现原文树拓扑的关键分支与 barcode gap 分布"（这两项对参数远不敏感），并在论文中明确写出 MOTU 数未能复现的事实——这个失败本身就是研究内容，它直接支持本课题"方法结果不稳定"的动机。
- 第 14 周末：完成 chronos vs BEAST2 对照。若 GMYC 的 MOTU 数在两种超度量化下相对偏差 > 25%，立即降级：把 GMYC 从四方法中剔除，主分析改为 ASAP / ABGD / mPTP 三方法，并把"GMYC 对超度量化方式高度敏感"作为一条独立结论报出。主结论框架不变。
- 第 26 周：完成主回归。若留一 R² < 0.15，按风险声明写"不可预测"结论，并把检验力估计与外部验证的失败率作为核心证据。
- 第 36 周结果冻结，第 44 周英文 PPT 初稿。
- 需提前核实而非边做边发现：(a) BOLD 批量下载的接口、配额与许可条款（能否在论文中再分发序列 ID 列表）；(b) 候选复现论文的补充材料是否给出可复原的登录号；(c) IQ-TREE 2 在 3,000 条序列上的内存峰值（须实测，16 GB 是否够）。三项在第 6 周前查清。
- 选择前提：适合愿意花时间打磨命令行管线、并接受"主结论可能是一条否定判据"的学生。这条路线的最大风险是数据清洗工作量（污染序列、NUMT、错误鉴定），须把清洗规则当成研究内容而非杂活来写。
- 预算裁剪顺序：先砍地理跨度属性（BOLD 坐标缺失率高），再砍稀释实验（H2），再砍外部验证集（从 10 个降到 5 个）。砍到只剩"30 个数据集 + 三方法 + 五个属性回归"时，主结论仍成立。