复现 Cao 人胎儿细胞图谱:四百万细胞之后,证据该回到哪一个统计单位?
一套包含四百多万个细胞的人胎儿图谱,很容易给人一种直觉:数据量如此庞大,组织差异和发育变化理应非常清楚,统计功效也不会成为问题。这个直觉把“细胞数量”与“生物学重复数”混在了一起。图谱里有 4,062,980 个细胞,承担人群层面重复的却是 28 个胎儿个体;在同一组织、同一细胞大类的配对比较中,实际可用的胎儿通常只有 4–7 个。细胞很多,可以让每个胎儿内部的表达均值估得更稳定,却不能凭空增加独立胎儿。
本文复现的重点由此落在一个很具体的问题上:完整读取数百万细胞的原始计数矩阵以后,哪些结论仍能在胎儿级统计单位、小样本精确检验和大规模多重校正下成立?围绕这个问题,分析依次检查全量计数闭合、胎儿级 pseudobulk 分母、跨组织配对、发育日龄关联、Tau 组织特异性、技术分区一致性,以及作者提供的 S4–S8 汇总资源。零发现、未定义值和不可识别的设计因素都保留在证据链里。
分析完整读取了 63,561 × 4,062,980 的 raw exonic-count loom,即 63,561 个带版本号的 Ensembl features 和 4,062,980 个细胞。矩阵中的精确计数总和为 4,295,529,686 counts,检测事件总和为 2,631,157,163。全量对象来自 28 个 fetuses、15 个 organs 和 77 个作者定义的 broad cell classes。本文所有群体推断都把 fetus 作为生物学重复;细胞、实验批次分区和可视化抽样均不充当独立重复。
这套分析得到的结果并不适合压缩成一张“显著基因榜”。组织配对在 family-level BH 校正后为 0 个 primary calls;发育日龄的精确全 family 校准也得到 0 个 primary discoveries。与此同时,大量方向在 leave-one-fetus-out 后保持稳定,Tau 排名显示连续的组织偏好,S8 中有少数方向出现同一图谱内的排名富集。把这些层次分开,才能同时保留有用线索与统计边界。
论文和数据来源
| 项目 | 来源 |
|---|---|
| 论文 | Cao et al. A human cell atlas of fetal gene expression,Science,2020 |
| DOI | 10.1126/science.aba7721 |
| PubMed | PMID 33184181 |
| GEO | GSE156793 |
| 主要计数输入 | S3 raw exonic-count loom |
| 细胞与基因元数据 | S1 cell metadata、S2 gene metadata |
| 作者汇总资源 | S4/S5 tissue expression/detection、S6/S7 organ-class expression/detection、S8 differential-expression resource |
原论文采用三级组合索引,对 15 个器官、超过 110 份组织样本进行表达图谱构建,最终获得约四百万个单细胞,并结合文献与其他图谱注释数百种细胞类型和亚型。论文重点包括广泛分布细胞类型的器官特化、胎儿红细胞生成部位以及与小鼠发育图谱的整合。这里的“超过 110 份样本”不能直接换算成 110 个独立胎儿;本次可核验的生物学个体为 28 个胎儿。
本文聚焦公开表达矩阵能够支持的胎儿级组织与日龄证据。作者的精细亚型体系、伴随的染色质可及性图谱、与小鼠图谱的整合,以及原文每一条器官生物学结论不在本文的完整复现范围内。GEO 提供的是外显子计数 loom;原始测序 reads 受控存放,因此本文也不声称完成 demultiplexing、比对和 read-level QC。
矩阵层、标签来源与连接键
S1 细胞表的 4,062,980 行与 loom 列顺序逐项一致,胎儿、器官、发育日龄、作者 broad class、assay 和 experiment batch 均无顺序错配。S2 基因表的 63,561 行也与 loom 行属性一致。带版本号的 Ensembl gene_id 在这 63,561 个 features 中唯一,因此所有矩阵连接和 S4–S8 对照都以它为键。
短基因名只用于显示。S2 中有 2,986 个 short symbols 对应多个带版本号的 Ensembl ID;直接按 symbol 合并会把不同 features 压成一行,尤其会污染 top-gene、marker overlap 和 held-out validation。文章后面提到的候选基因集合都遵守同一连接规则,不用短名称承担唯一性。
标签来自作者元数据中的 broad cell classes。本文没有用最高分规则再造一套细胞身份,也没有把 broad class 当成现代细粒度 taxonomy 的替代品。公开 loom 的 Size_Factor 字段对全部细胞缺失,因而不能宣称复用了作者的 size-factor normalization。定量层使用精确计数构建 CPM 和检测比例;组织效应在胎儿内比较 log2(CPM+1),Tau 使用各器官的等胎儿权重均值。
证据路线:从工程闭合走到胎儿级推断
分析路线可以压缩为六个相互约束的层次:
- 扫描完整 loom,逐 block 同时累加 count 与 detected-cell events,并核对原矩阵、分组矩阵、细胞总计和 H5 汇总层是否闭合。
- 将 1,378 个技术单元合并为 1,178 个胎儿×器官×日龄×作者 broad class 的 biological pseudobulks。
- 预先固定描述阈值、推断阈值、共享 class、组织配对和日龄 stratum 的资格规则。
- 对同一胎儿内的组织配对做精确检验、family-level BH 和 leave-one-fetus-out 方向检查。
- 对日龄关联做全 family 的精确 day-label permutation;对组织特异性另行计算无显著性含义的连续 Tau。
- 用技术分区和作者 S4–S8 资源检查一致性,同时维持 assay、同图谱 held-out 和可视化抽样的边界。
以下 12 张图各自只回答一个问题。每张图后都同时说明它支持的证据和不能由它推出的结论。
从 1,378 个技术组到 1,178 个生物学组
全量扫描先按照胎儿、器官、发育日龄、作者细胞大类、assay 和 experiment batch 六个字段建立 1,378 个技术组。这个立方体保留了数据如何被实验分区的信息,适合检查每个分区的计数、检测比例和细胞数,却不适合直接进入群体统计。若同一胎儿、同一器官、同一日龄和同一细胞大类被拆到两个实验批次,将两列当成两个样本会人为增加自由度。
生物学聚合因此去掉 assay 与 experiment batch 两个技术维度,把属于同一 fetus×organ×day×class 的分区相加,得到 1,178 个 biological pseudobulks。计数以无符号整数求和;检测分子使用各分区中表达大于零的细胞数求和,检测分母使用各分区总细胞数求和。这里没有先算每个分区的检测比例再做简单平均,因为大小不同的分区若获得同等权重,会改变真实的细胞级检测分母。
聚合前后计数总和都为 4,295,529,686,检测事件总和都为 2,631,157,163,说明技术层到生物学层没有丢失或重复。每个 pseudobulk 仍然只是一个胎儿在某一器官和细胞大类中的观测;同一胎儿贡献的多个 pseudobulks 在统计上有关联。组织分析通过胎儿内配对利用这种关联,日龄分析则明确承认胎儿与日龄无法拆分。
图1要问:完整 raw loom 是否真的被全量读取并精确闭合?

图1|raw exonic-count layer 的全量扫描与精确闭合。 左侧列出 63,561 个 features、4,062,980 个 cells、1,378 个技术分组、994 个连续 feature blocks、4,295,529,686 个 count events 和 2,631,157,163 个 detection events;右侧显示 count 与 detection 的权威汇总减独立总计均为 0。闭合证明完整计算,不构成生物学发现。
全量扫描逐 block 读取所有 4,062,980 个细胞,没有用少量 marker streaming 或抽样矩阵代替定量分析。994 个 block 中有 88 个全零 feature blocks,它们具有相同压缩内容并不表示 block 被重复或漏读;连续起止位置、文件尺寸、矩阵形状、非负整数、原始计数到分组计数、细胞检测事件到分组检测事件均通过闭合检查。扫描自然完成,wall time 为 2:56:48,peak RSS 为 9,799,756 KiB,退出状态为 0。这些运行信息只用于证明矩阵确实被处理。
权威计数层采用 uint64 精确聚合。float32 分组聚合在数值检查中曾于 8 个 blocks 出现计数偏差,最大偏差为 70,905 counts,因此该数值层不参与任何结果解释。最终 count、library size、逐细胞总计和来源汇总都精确等于 4,295,529,686;detection 两条独立求和路径都等于 2,631,157,163。
图1支持的是工程闭合与矩阵层边界。它不能说明某个器官更特殊,也不能把四百多万个细胞换算成四百多万个独立重复。要进入生物学推断,还要先把技术分区还原到胎儿级单位。
图12要问:多少胎儿级 pseudobulks 达到描述和推断阈值?

图12|胎儿级分母与资格阈值。 1,178 个 biological groups 定义为 fetus×organ×development day×main class,其中 874 个达到 20-cell 描述阈值、730 个达到 50-cell gene-level inference 阈值。groups 是同一批 28 个胎儿的分层观测,不能当成 1,178 个独立重复。
1,178 组中,874 组达到 20-cell 描述阈值,730 组达到 50-cell 推断阈值。组内细胞数跨度极大:最小为 1,中位数为 109.5,第 25 百分位约为 18.25,最大达到 467,066。因而“全图谱细胞很多”与“每个胎儿—器官—细胞大类组合都有足够细胞”是两件不同的事。低于阈值的组继续保留在分母表里,标为不可估计或仅作描述,没有被补成零效应。
这里的 1,178 组也不能当成 1,178 个独立样本。同一个胎儿会贡献多个器官和多个 broad classes,独立生物学个体仍为 28。20 与 50 是预先固定的组内细胞门槛,用来控制均值和基因层估计的可用性;它们不改变胎儿是重复的事实,也不构成结果显著性。
为什么不可估计的组合不能填成零?
1,178 个生物学组中有 304 组低于 20 个细胞,另有 144 组处于 20–49 个细胞之间。前 304 组连稳定的描述性均值都缺少足够支持,后 144 组可以进入器官/class 分布等描述,却不进入基因层推断。若把这些组的效应、表达或检测结果统一填成零,会把“数据不足”改写成“测得无效应”,同时人为扩大分母。
保留不可估计状态对组织配对尤其重要。一个胎儿在 organ A 有充足细胞、在 organ B 没有同一 class,并不构成 organ B 的零表达配对;这个胎儿应从该 pair 的基因检验中退出,同时继续出现在资格与缺失分母表里。日龄 stratum 也采用相同原则:未达到胎儿数、日龄数或 50-cell 门槛的组合保留为描述性记录,不能补入零斜率。
这种处理会让可检验样本数变小,却保留了观测与缺失的区别。后续实验设计可以据此识别哪些器官、class 或日龄缺少胎儿覆盖,而不会被一大片人工零值误导。图2展示的共享资格正是建立在这条分母纪律之上。
图2要问:哪些 broad classes 具备跨器官、多胎儿的支持?

图2|共享 broad class 的资格图。 每个格子的数字表示在该器官和 broad class 中达到 20-cell 描述阈值的胎儿数。共享 class 需要总计至少 5 个胎儿、至少 3 个合格器官,并且每个纳入器官至少有 3 个胎儿支持。格子表示 fetus support,不表示 abundance、细胞比例或独立重复数。
77 个作者 broad classes 中有 11 个通过共享规则。支持面最广的是 vascular endothelial cells,覆盖 13 个合格器官;stromal cells 覆盖 11 个;myeloid 与 lymphoid cells 各覆盖 10 个。astrocytes、lymphatic endothelial cells、mesothelial cells 和 Schwann cells 各有 3 个合格器官,刚好达到跨器官分析的最低结构要求。热图同时显示,某些 class 在个别器官中细胞很多,却未必在多个胎儿和多个器官中均衡出现。
这张图用于筛选后续可比较的 class 和 organ 组合。颜色与数字表达胎儿支持,不表达总体细胞比例、胎儿人群中的发生率或器官丰度。一个格子里有 11 个胎儿,也不能与另一格子的 3 个胎儿简单比较为“表达更强”;这里没有展示表达量。
图3要问:4–7 个配对胎儿能达到怎样的精确 P 值分辨率?

图3|小样本精确检验的离散分辨率。 44 个 tissue pairs 中有 17 对由 4 个胎儿支持、15 对由 5 个胎儿支持、7 对由 6 个胎儿支持、5 对由 7 个胎儿支持。在无零差值的理论边界下,多重校正前的最小双侧 P 依次为 0.125、0.0625、0.03125 和 0.015625;虚线为 0.05。功效有限不能证明组织没有差异。
组织比较使用同一胎儿内、同一 broad class 的配对 log2(CPM+1) 差值。配对设计同时控制胎儿身份和该胎儿对应的日龄,比跨胎儿直接比较器官均值更稳妥。代价是可参与某个器官对的胎儿只剩 4–7 个。n=4 或 n=5 时,理论上连未校正的双侧 P<0.05 都无法达到;n=6 或 n=7 虽能跨过 0.05,之后仍要面对一个由 63,561-feature universe 衍生、每对实际含 18,622–35,047 个可估计基因的多重检验 family。
图3展示的是数学可达范围,不对应观测到的显著性,也没有证明组织效应不存在。它解释了为何“所有胎儿方向一致”仍可能无法形成校正后发现,并为图4中方向稳定与 inferential call 的分离提供背景。
图4要问:方向稳定的行很多,为何 primary tissue calls 仍为零?

图4|方向稳定与多重校正是两个证据层。 44 个 tissue pairs 各有 4–7 个配对胎儿;左图给出 leave-one-fetus-out 后均保持同一均值方向的基因数,右图给出 class×organ-pair family 内 BH q≤0.05 的调用数,结果为 0。零调用不能解释成零组织差异。
44 个配对共有 1,201,791 个可估计 tested rows,其中 680,990 行在完整数据和每一次可估计的 leave-one-fetus-out 分析中保持相同方向。单个 pair 的方向稳定行介于 8,145 与 24,790 之间;例如 stromal cells 的 Heart–Lung、Lung–Muscle 和 Adrenal–Lung 比较分别出现 24,790、24,040 和 22,807 个方向稳定行。如此大的数量说明,稳定方向本身是很宽松的描述条件。
同一批数据在 family-level BH 后,44 个 pairs 的 q≤0.05 基因数全部为 0,stable_tissue_specific_call 也为 0。这里没有用未校正 P、效应大小排序或 LOO 稳定性替代 primary gate。方向稳定的行可以用于提出后续验证候选,却不能写成组织差异发现。
0 个 primary calls 也不等于胎儿组织之间没有差异。它表示在 4–7 个配对胎儿、离散精确 P 和每个 family 大量 features 的条件下,没有基因跨过预先规定的多重校正门槛。扩大独立胎儿数、预先限制假设家族,或在独立队列中验证排序,才可能提高可解释的推断功效。
组织配对结果需要拆成五个字段阅读
每个基因的组织比较至少包含五类信息:配对胎儿数、完整数据的均值与中位数差、正向/负向/并列胎儿数、精确 P 与 family-level BH q,以及 leave-one-fetus-out 后的效应范围。效应大小回答“差多少”,方向计数回答“有多少胎儿同向”,精确 P 回答“在配对零假设下这种符号与秩组合多罕见”,BH q 处理同一 class×organ-pair family 中的大量基因,LOO 则检查结论是否由某一个胎儿支撑。
这些字段不能相互替代。均值差很大可能来自一个极端胎儿;方向一致可能因为 n 很小;LOO 同向可能与未校正 P 同时出现,却仍无法通过数万行的 BH。反过来,若将完整矩阵先按效应挑出少量基因,再只对这批候选报告 P,检验 family 会被事后缩小,显著性也会偏乐观。本分析保留每个 pair 的完整可估计集合,0 个 BH calls 因而具有明确分母。
1,201,791 个 tested rows 是 44 个 families 的行数之和,同一个 Ensembl feature 可以在多个 class 或器官对中重复出现。680,990 个稳定行也按 pair-specific rows 计数,不能理解为 680,990 个唯一基因。公开候选表若要进一步压缩,应同时保留 pair、class、胎儿数、检测门槛和原始 family 信息,避免把跨 pair 重复当作独立证据。
图5要问:跨胎儿日龄关联能否建立发育连续性?

图5|发育日龄的小样本精确校准。 左图比较探索性渐近法与 primary exact day-label permutation 在每个 stratum 中得到的 BH q≤0.05 行数;右图展示各 stratum 的最小精确置换 BH q。精确 primary discoveries 为 0,渐近法产生的 223 行只承担校准审计角色。每个 fetus 只有一个 day,fetus 与 day 不可分别识别。
日龄分析包含 51 个 organ×main-class strata、1,409,275 个 tested rows。每个 stratum 至少有 6 个胎儿、4 个不同日龄,并要求每个胎儿组至少 50 个细胞;实际胎儿数分布为 6–11。每个 stratum 内对 day labels 做双侧精确置换,再在该 stratum 的完整 tested family 中进行 BH。全体 exact BH<0.05 行数为 0,联合 primary call 也为 0。最小原始精确置换 P 可达到约 3.21×10^-6,但全 family 校正后的最小 q 仍约为 0.1235。
渐近 Spearman tail 在这些小 n strata 中给出 223 个 BH<0.05 行;若继续叠加方向一致、leave-one-fetus-out 和固定日龄分箱等条件,只会在同一探索性链中进行复合筛选,不能产生独立证据。精确全 family 校准没有支持任何一行进入 primary 发现集,因此这 223 行不能作为 stage gene list。它们的用途是展示小样本渐近尾部可能过于乐观。
更根本的设计限制来自胎儿与日龄的绑定:每个 fetus 只有一个 development day。数据中没有同一胎儿跨日龄的重复,也没有足够的独立结构把个体差异与日龄差异完全拆开。51 组结果只能称为 cross-fetus stage association。固定在全局日龄三分位点的 early/middle/late 分箱仅作敏感性检查;UMAP 排列、分箱均值或相关系数都不能建立轨迹、谱系方向和发育机制。精确 stage 主发现为 0,也不等于人胎儿组织没有随日龄变化。
为什么精确日龄置换必须覆盖完整 tested family?
日龄关联中,每个基因先在一个器官×细胞大类 stratum 内形成胎儿级表达向量,再把观测到的 day labels 在这些胎儿之间置换。双侧精确 P 来自所有允许的标签排列,不依赖大样本 t 分布近似。不同基因可能共享同一缺失模式和秩结构;全体 51 strata 中归纳出 364 种独特的 gene-rank null patterns,可以在维持精确性的同时复用计算。
置换范围与 BH family 都在查看结果前固定。若只对渐近 P 最小的 223 行补做 permutation,再对这 223 行校正,就会忽略它们从 1,409,275 个 tested rows 中被筛出的过程。这里对每个 stratum 的全部可检验基因产生精确 P,再在该 stratum 内完整执行 BH。即使某个基因的原始精确 P 很小,只要同 family 中同时检验的行数很大,q 仍可能高于 0.05。
图5右侧显示的最小 q 大多远离 0.05,说明零发现来自精确尾部与多重校正的共同结果,不能归因于机械设置了过严阈值。若未来只关心少数预注册通路,可以在新数据采集前定义更小且有生物学依据的检验 family;不能在当前数据中根据结果反向挑选 family,再把探索性结果改名为主发现。
图6要问:Tau 如何描述共享细胞大类的组织偏好?

图6|Tau 是等胎儿权重的连续描述量。 对每个 broad class,先在每个器官内对胎儿等权计算 CPM 均值,再定义 Tau = Σ(1 - organ_mean/max_organ_mean)/(number_of_organs - 1)。图中的线段和箱体表示 Tau 的 5–95%与25–75%区间,红点为中位数。Tau 没有 P 值或 inferential null。
Tau 接近 0 表示在合格器官间较均匀,接近 1 表示表达集中于较少器官。11 个 classes 的中位 Tau 从 lymphatic endothelial cells 的约 0.623 到 erythroblasts 的约 0.925;vascular endothelial、myeloid、lymphoid 和 stromal cells 的中位数约为 0.83–0.85。这样的分布提示,同一个作者 broad class 跨器官出现,并不意味着其表达程序在各器官相同。
Tau 没有预先登记的 inferential null,也没有 P 值。高 Tau 不能自动转换成“显著器官特异基因”,更不能从中推出功能或机制。它适合排列候选、比较 class 间分布,并与 detection、胎儿支持和 leave-one-fetus-out 稳定性一起阅读。其分母问题由图7单独处理。
图7要问:Tau 的稳定比例为何必须声明分母?

图7|Tau 的 defined 与 all-feature 两套分母。 primary denominator 为 335,595 个 defined-Tau rows;699,171 个 all-feature rows 只作敏感性分母,其中 363,576 个 all-zero rows 的 top/second organ 标签保存为 NA。左图报告 LODO 与 expression/detection 一致性,右图展示分母切换造成的位移。
11 个 classes 合计有 699,171 个 gene-class rows。335,595 行在至少一个合格器官中具有非零表达,因此 Tau 有定义;另有 363,576 行在所有合格器官中的 expression 和 detection 都为零。对全零向量,Tau 公式的分母 max_organ_mean 为 0,top organ 也没有生物学含义。所有这类行的 top/second expression 与 detection 标签均为空,没有让数组 argmax 随意指定第一个器官。
在 335,595 个 defined-Tau rows 中,各 class 的 LODO top-organ 稳定比例为 21.7%–63.5%,中位数约 38.4%;expression 与 detection 的 top-organ 一致比例为 62.9%–91.2%,中位数约 81.0%。若把 363,576 个全零行放回 699,171 的全 feature 分母,两类比例都会明显下降。独立数值复核显示,分母切换造成的最大绝对位移可达约 0.329 和 0.626。
图7支持的是“候选排名对遗漏单个胎儿有多稳”“表达峰值与检测峰值是否一致”以及“分母选择会改变多少”。LODO 稳定不提供显著性检验,defined-Tau 比例也不能称为全 atlas 无条件比例。排名较稳的 features 可以进入独立实验候选清单,但在独立样本验证前不应称为器官 marker。
为什么表达强度与检测比例要并行报告?
CPM 与检测比例描述表达的两个侧面。某个基因可能只在少数细胞中表达很高,形成较大的 pooled CPM,却在大多数细胞中检测不到;也可能在许多细胞中稳定低表达,检测比例较高而总计数不突出。单独使用其中一个指标,容易把稀疏高表达与广泛低表达混为同一类组织偏好。
本分析在每个 fetus-organ-class 组内保存精确 count、library size、positive-cell numerator 和 cell denominator,再对器官均值采用胎儿等权。这样,大胎儿样本中的大量细胞不会让该胎儿在 Tau 中获得更高权重。top-expression organ 与 top-detection organ 一致时,只能说明两种描述在器官排序上同向;一致比例较高仍不提供显著性,也不能排除细胞组成或检测深度造成的影响。
图7中 expression/detection 一致比例在 defined-Tau 分母下多为 0.63–0.91,高于 LODO top-organ 稳定比例。这个差别提示,表达与检测可以在完整数据中给出相同器官排序,但遗漏一个胎儿后 top organ 仍可能变化。候选优先级应同时考虑两项指标:表达/检测一致可减少度量选择带来的歧义,LODO 则暴露个体依赖性。两项都通过的候选仍需独立样本确认。
图8要问:技术分区一致性会不会被低细胞数拖累?

图8|technical partitions 是同一生物学组内的技术描述。 159 个 multi-partition groups 中,136 个每个 partition 均≥20 cells,23 个含小 partition;partition 不充当生物学重复。横轴为最小分区细胞数,纵轴和颜色分别表示表达与检测一致性。0 个组含配对双 assay,因此 assay sensitivity 不可识别。
共有 159 个 biological groups 含多个 experiment-batch partitions。136 组的每个 partition 都至少有 20 个细胞,23 组至少包含一个小于 20 个细胞的 partition。前一组的最小 partition 细胞数中位数为 666.5,pairwise expression minimum Spearman 中位数约 0.857;小 partition 组的最小细胞数中位数只有 4,相同相关指标中位数约 0.363。极端情况下,一个只有 1 个细胞的 partition 使表达与检测相关都降到约 0.027。
这张图说明低细胞技术分区会显著削弱一致性,因此只报告全体中位数会掩盖风险。159 组中的 technical partitions 仍属于同一个胎儿×器官×日龄×class 单元,它们不能增加生物学样本量,也不能用于细胞级伪重复。
assay 层还有更严格的限制:159 个多分区组全部只有一种 assay,整个数据中 0 个 biological groups 同时包含两种 assays。因此 within-group assay sensitivity 不可识别。这里不能写“不同 assay 结果一致”,也不能把缺少配对设计解释成 assay effect 为零。
图9要问:全量 loom 汇总与作者 S4–S7 资源一致到什么程度?

图9|detection 可精确重构,expression 只作 estimator-boundary 下的秩对照。 187 个 S4–S7 columns 的 detection 在 1×10^-12 容差内全部精确重构;expression 与 pooled exact-count CPM 的 Spearman 最低约 0.99684。高秩一致不表示两个 expression estimators 数值相同。
复核覆盖 15 个 S4/S5 tissue columns 和 172 个 S6/S7 organ-main-class columns,共 187 列。S5/S7 的 detection fractions 在全部列中都能从 full loom 精确重构到 1×10^-12 容差内,最大绝对差为 1.11×10^-16。这个结果为细胞顺序、器官/class 聚合和正表达计数提供了强一致性证据。
S4/S6 的 expression summaries 与 pooled exact-count CPM 的 Spearman 范围为 0.9968376–0.9999529,排名高度一致。最低相关出现在只有 36 个重构细胞的 Cerebrum–Megakaryocytes 列,但相关仍约为 0.99684。高相关不代表两个 estimator 数值相等;作者表达汇总与 pooled CPM 的计算细节可能不同,本文只报告 rank agreement 和重叠层面的复核。
图9也不能替代前面的胎儿级组织推断。作者汇总列通常把多个胎儿合在一起,适合检查资源一致性,却不保留每个组织配对所需的独立胎儿结构。检测精确重构增强了数据链可信度,不能把图4的 0 primary tissue calls 改写成显著结果。
图10要问:secondary tissue 排名是否在 S8 中超过大小匹配的 null?

图10|S8 只能作为同一 atlas 内的 held-out 排名校准。 组织效应的 primary gene set 在查看 S8 之前已经固定,S8 没有参与 tissue-effect construction、filtering 或 tuning。分析把 44 个 tissue pairs 展开为两个 favored directions,共 88 个 pair-directions;每个方向按预先规定的 expression/detection、LOO 和效应排序规则取不超过 100 个 candidates,再进行 10,000 次大小匹配的 null 抽样。横轴是 null 下平均重叠数,纵轴是实际 S8-supported genes;颜色是 raw one-sided permutation P,跨 88 个方向校正后的 q 需另行计算,颜色也不表示效应大小或验证概率。
88 个方向中有 4 个 nominal P<0.05,其中 2 个在 across-88 BH 后 q<0.05;69/88 个方向的 candidate universe≤top-N,属于退化或信息受限方向。primary tissue-effect calls 仍为 0。这四组数字必须一起阅读,省略其中任何一组都会把 secondary calibration 写成 discovery。
两个 across-88 BH 方向都来自 vascular endothelial cells 的 Cerebrum-favored 排名:Cerebrum–Heart 方向实际重叠 80,大小匹配 null 均值约 64.50,raw P≈1.0×10^-4、BH q≈0.0088;Cerebrum–Muscle 方向实际重叠 77,null 均值约 65.89,raw P≈0.0010、BH q≈0.0440。另外两个 nominal 方向为 Cerebellum–Muscle 中 Cerebellum favored,以及 Cerebrum–Liver 中 Cerebrum favored,跨 88 方向校正后均未达到 0.05。
这些点提供的是同一 atlas 内的 ranking clue。S8 与 primary 计数来自同一数据资源,缺少独立胎儿队列和正交实验;69 个退化方向又限制了 top-N null 的信息量。它不能称为外部验证、组织方向验证,也不能挽救、替代或扩大 0 个 primary tissue calls。连接仍以带版本号 Ensembl ID 完成,short symbol 只用于显示。
S8 的 secondary candidates 应该怎样使用?
进入 S8 校准的候选遵循预先固定的方向规则:完整数据和所有可估计 LOO 均保持同一效应符号,expression 与 detection 指向同一 favored organ,且 favored organ 的平均检测比例至少为 25%。满足条件的基因按配对平均 log2(CPM+1) 差值绝对值排序,每个方向最多取 100 个,带版本号 Ensembl ID 用作确定性并列键。88 个方向合计产生 6,342 个方向特异的 candidate slots;这个数字包含同一基因在不同 pair 或方向中的重复,不代表 6,342 个唯一发现。
大小匹配 null 从同一方向的完整 candidate universe 中反复抽取相同数量的基因,回答“排序前端的 S8 support 是否多于随机挑取同规模候选”。当 candidate universe 本身不超过 top-N 时,抽样集合与全集高度重合,观测重叠会接近 null 重叠,检验缺少区分能力。69 个退化方向因此需要单独计数,不能归为 69 个阴性验证。
两个跨 88 方向 BH 的点可用于安排外部复核优先级,例如优先查看脑组织 vascular endothelial cells 的候选是否在独立胎儿、原位表达或蛋白层重复出现。合理用途是“决定先验证谁”,不包括“用同一图谱证明组织方向成立”。若新的独立数据改变了候选 universe、细胞定义或检测深度,还需要重新计算排序和 null,不能直接搬用这里的 P 与 q。
图11要问:作者全局 UMAP 能提供什么背景?

图11|rendering-only 的全局组织背景。 图中展示 28,124 个确定性抽样细胞,抽样覆盖全部 1,178 个 fetus-organ-class strata,每个 stratum 上限为 30 个细胞;坐标直接来自作者 S1 global UMAP。抽样不参与推断,也不替代 4,062,980 个细胞的完整定量分母。
器官颜色在部分区域形成明显块状结构,在另一些区域交叠。这个画面帮助读者理解图谱同时包含器官特异群体和跨器官共享群体,也解释了为何后续要在 broad class 内做器官比较。所有数值结果仍来自 4,062,980 个细胞的完整聚合或胎儿级表,抽样文件没有进入 tissue、stage、Tau、technical 或 S8 推断脚本。
UMAP 只能提供背景。点的接近、分支形状、颜色梯度和局部混合不能推出发育轨迹、细胞命运、细胞相互作用或机制;图中抽样比例也不能用于估计器官或 class 丰度。若要研究连续发育,需要独立的轨迹模型、足够的跨日龄胎儿和对个体效应的明确处理。
与原论文和作者资源的对应关系
原论文将该数据定位为人胎儿多器官表达图谱,并强调广泛分布细胞类型的器官特化。本文的 11 个共享 broad classes、44 个同胎儿组织 pairs 和 Tau 分布与这个问题直接相连,但采用了更严格的统计拆分:器官特化线索可以出现在效应方向、Tau 或 S8 排名中,群体发现则必须回到胎儿级配对与多重校正。
三类对应关系需要分开:
- 可直接复核的资源事实:S3 loom 的 63,561×4,062,980 形状、S1/S2 顺序、28 个胎儿、15 个器官、作者 broad class、S5/S7 detection fractions,以及作者 global UMAP 坐标。
- 本文增加的统计层:1,178 个 biological pseudobulks、20/50-cell 资格门、44 个同胎儿 tissue pairs、51 个精确 stage families、等胎儿权重 Tau、technical partition size sensitivity 和 S8 大小匹配 held-out ranking calibration。
- 没有建立等价关系的部分:作者全部精细细胞亚型、染色质可及性图谱、原论文所有 differential gene lists、与小鼠发育图谱的整合、read-level 流程,以及作者 expression estimator 与 pooled exact-count CPM 的数值同一性。
原论文摘要还指出胎儿红细胞生成可见于肾上腺等部位。图2确实显示 erythroblasts 在多个器官具有胎儿支持,图3与图4也包含 Adrenal–Liver 的 erythroblast 配对;但本文没有以原论文该命题为预注册的独立检验,也没有用正交组织学确认,因而不把这一对应写成新的机制复现。
如果把细胞量误当成样本量,会出现哪些具体偏差?
以 organ A 和 organ B 的一个共享 class 为例,假设四个胎儿各自贡献数万细胞。若把每个细胞作为独立观测,检验会把同一胎儿内共享的遗传背景、日龄、取材条件和实验过程当成大量独立信息,标准误会被压得很小,几乎任何微小差异都可能得到醒目的 P 值。细胞级图看起来很稳定,却没有回答差异能否跨胎儿重复。
胎儿级 pseudobulk 将每个胎儿在每个器官中的细胞整合为一个表达与检测观测,再做同胎儿配对。数万细胞的作用是降低该胎儿均值的测量噪声,检验自由度仍由四个胎儿决定。技术分区同理:把同一胎儿的两个 experiment batches 当成两个胎儿,会同时放大样本量并低估批次相关性。
这种伪重复还会影响候选验证。若组织差异基因由细胞级 P 选出,再用同一批细胞的作者汇总表或 S8 检查 overlap,构建与验证共享了相同个体和大量相同技术噪声。本文把 primary 选择放在胎儿配对层,并把 S8 限定为同图谱 held-out ranking calibration,仍然不能获得外部独立性,但至少没有把细胞数包装成胎儿数。
结论审计:哪些可以写,哪些只能保留为线索?
| 证据等级 | 可以保留的表述 | 需要拒绝的扩张 |
|---|---|---|
| 支持 | 完整 raw exonic-count loom 已全量读取,count/detection 精确闭合 | 细胞数量等同于独立重复数 |
| 支持 | fetus 是生物学重复;1,178 个 pseudobulks 中 874/730 个达到描述/推断阈值 | 把 pseudobulk、cells 或 partitions 当成 1,178 个独立个体 |
| 支持 | 44 个 tissue pairs、每对 4–7 个胎儿,family-level BH calls=0 | “胎儿组织之间没有差异” |
| 支持 | 51 个 stage families、1,409,275 tested rows,exact primary=0、joint primary=0 | 把 223 个渐近行写成 stage discoveries |
| 支持 | Tau 是等胎儿权重的连续描述量;335,595 defined rows 与 363,576 all-zero rows 必须分开 | Tau 显著、Tau 是二元发现或全零行有 top organ |
| 支持 | 159 个 multi-partition groups 中 136 个分区充分、23 个含小分区 | technical partition 是生物学重复 |
| 支持 | 0 个组同时含两种 assays,assay sensitivity 不可识别 | 不同 assays 已证明一致 |
| 支持 | 187 个作者列的 detection 精确重构,expression 高度秩一致 | 两种 expression estimators 数值相同 |
| 探索 | LOO 稳定的 tissue rows、稳定 Tau 排名可用于独立验证排序 | 把方向稳定写成校正后组织发现 |
| 探索 | S8 有 4 nominal、2 across-88 BH 的 ranking enrichment clues | S8 是外部验证或能恢复 primary tissue calls |
| 背景 | 28,124-cell UMAP rendering 展示全局器官结构 | 从 UMAP 推断轨迹、相互作用、命运或机制 |
深度解读:这篇复现应该怎么读?
这篇复现回答了什么?
它回答的是一个证据尺度问题:当完整四百万细胞对象被实际读取后,结论回到 28 个胎儿及其配对结构,会剩下哪些可审计的结果。答案包含两部分。工程与资源层非常牢固:计数闭合、检测闭合、元数据 crosswalk、作者 detection summaries 和全局嵌入都能复核。群体推断层受到小样本和设计限制:tissue 与 stage 的 primary calls 均为 0,但方向稳定、Tau 和 held-out ranking 仍提供后续候选。
raw matrix 与 processed resource 的边界在哪里?
主要计数输入是 raw exonic-count loom,属于未标准化的计数矩阵,与原始 FASTQ 的 read-level 数据层不同。组织与日龄分析从这套计数构建胎儿级 CPM、检测比例和 pseudobulk。S4–S7 是作者生成的 processed summaries,用来复核 detection 与 expression ranking;S8 是同一 atlas 内预留的 differential-expression resource,用来做 secondary ranking calibration。任何 processed resource 都没有替代 full loom 的分母和胎儿级推断。
图链应该按什么顺序阅读?
先读图1和图12:前者证明完整输入与数值闭合,后者确定胎儿级单位和资格阈值。图2说明哪些 broad classes 有跨器官支持。图3与图4合起来解释 tissue 结果,图3给出可达 P 分辨率,图4分开方向稳定与 primary calls。图5处理 stage 的渐近法偏差和 fetus/day 不可识别。图6、图7再读连续的 Tau 及其分母。图8检查技术分区,图9检查作者资源,图10只读 same-atlas ranking calibration。图11放在末尾当作地图背景,避免让 UMAP 形状主导推断。
哪些地方容易被过度解释?
最常见的误读是把“零发现”写成“零生物学效应”。tissue primary=0 受到 4–7 个胎儿、离散 P 和大 family 校正限制;stage primary=0 还叠加 fetus/day 绑定。另一个误读是把稳定性当显著性:LODO 只说明结果对遗漏一个胎儿的方向是否敏感,不能生成 P 值。Tau 的高值说明相对集中,不提供 inferential null。S8 与 primary 数据同源,无法承担外部验证。technical partitions 和可视化 cells 都不能替代胎儿。
如果继续深入,下一步做什么?
组织效应需要增加每个 pair 的独立胎儿数,并在采样前限定少量组织、class 和 gene-program hypotheses,避免把数万 features 同时交给极小 n。发育日龄研究需要同一日龄的多个胎儿、相邻日龄更均衡的覆盖,以及能把个体差异与日龄效应分开的设计。Tau 候选应进入独立胎儿队列、原位实验或蛋白层验证。assay 影响需要在同一 fetus-organ-day-class 中同时包含两种 assays 的桥接样本。S8 排名应在独立 atlas 中重新定义 support 集合后再检验。
这类超大 atlas 项目最容易在哪些环节出错?
把下载完成当成分析完成。 数十 GB 文件存在于磁盘并不能证明矩阵被读过。需要保留输入形状、连续 block 范围、原始到汇总的 count/detection closure、运行时间与内存记录。图1中的两个零 delta 比“文件已下载”更能说明全量执行。
把 cells、libraries、partitions 和 fetuses 混为一谈。 四百万细胞提升个体内均值精度,28 个胎儿决定群体推断上限。1,378 个技术组压缩成 1,178 个 biological groups 后,仍不能说有 1,178 个独立重复。统计单位必须由实验设计定义。
在小 n 中默认渐近 P 足够准确。 图5显示渐近法可给出 223 个 BH 行,而全 family 精确标签置换为 0。只在一小部分候选上补做 permutation 也不够,因为 BH 的分母必须覆盖预先规定的完整 tested family。
让未定义值偷偷进入分母。 363,576 个 all-zero gene-class rows 没有 Tau,也没有 top organ。若用数组位置替代 NA,会制造大量看似明确的器官标签;若将它们混入主分母,又会显著压低稳定比例。图7把两套分母并排展示,避免这类错误隐身。
把同源资源称为外部验证。 S8 没参与 primary tissue construction,因此具有 held-out 意义;它仍来自同一 atlas,不能消除样本来源、批次和胎儿结构的共享。held-out 与 external 是两种不同的独立性。
可复用的方法学经验
这套证据链适用于许多大型公共单细胞图谱。一个可审计的流程至少要保留以下要素:
- 明确矩阵层是 raw counts、processed counts 还是作者汇总表,并区分 count matrix 与 raw sequencing reads。
- 以唯一 feature ID 连接矩阵和外部资源,基因 symbol 只负责显示。
- 对全量矩阵同时保存 count 与 detection closure,抽样只用于渲染。
- 在任何检验前确定真实生物学重复、描述阈值、推断阈值和完整多重检验 family。
- 对小样本优先检查精确 P 的可达分辨率,并让 permutation 覆盖整个 tested family。
- 将效应大小、LODO 稳定性、显著性和 held-out ranking 分成不同字段,避免用其中一层替代另一层。
- 对连续 specificity 指标公开公式、未定义条件、主分母与敏感性分母。
- 对 technical partitions 只做组内一致性描述;没有 paired assay 时明确写不可识别。
- 每张图保留来源表、生成代码、参数和输入哈希,使图上的每个数字都能回到计算表。
这些要求会让文章少一些醒目的“发现数”,却能让候选排序、阴性结果和设计缺口都具有下一步实验意义。对于已有大型 atlas 或积压多年的单细胞项目,先建立这样的分母与可追溯链,通常比继续增加 UMAP、marker 列表和细胞级 P 值更有效。
局限性
本文保留全部 4,062,980 个公开细胞,但分析起点仍是 GEO 释放的 exonic-count loom,无法覆盖受控 raw reads 的 demultiplexing、比对与 read-level QC。作者 broad classes 直接来自元数据,本文没有重新建立数百个精细亚型,也没有复现伴随的 chromatin accessibility atlas。
组织配对只有 4–7 个胎儿,stage strata 只有 6–11 个胎儿。每个胎儿只对应一个 development day,使个体与日龄不可分别识别。Tau 没有 inferential null;S8 多数方向退化且与 primary 分析同源;technical assay 对比缺少同组双 assay。上述限制决定了候选排名可以保留,机制、轨迹、外部验证和“无差异”结论都不能从现有证据推出。
小结
Cao 人胎儿图谱的公开表达对象确实能够支持一次完整的四百万细胞全量计数复现:63,561 个 features、4,062,980 个 cells、4,295,529,686 counts 和 2,631,157,163 个 detection events 全部闭合。进入生物学推断后,重复单位回到 28 个 fetuses。44 个同胎儿组织 pairs 在 BH 后有 0 个 primary calls;51 个 stage families 在精确全 family 置换后也有 0 个 primary discoveries。Tau 提供连续的组织偏好描述,但必须拆分 335,595 个 defined rows 与 363,576 个 all-zero rows。S8 的 4 个 nominal、2 个 across-88 BH 和 69 个退化方向只构成同一 atlas 内的排名校准,不能改变 primary tissue=0。
这篇复现留下的主要成果是一条完整、可核验且允许失败的证据链:全量矩阵经过了真实计算;胎儿级分母没有被细胞数量稀释;零发现没有被包装成无效应;描述性排名也没有被提升为机制。类似的大型图谱若要形成可进入论文或报告的结论,可以沿同样顺序逐层收紧问题、统计单位和验证边界。