复现 Liao 新冠 BALF 单细胞数据：细胞很多，证据仍要回到 13 位供者

Liao 等人在 2020 年发表的新冠肺炎 BALF 单细胞研究，记录了早期 COVID-19 队列的肺泡灌洗液免疫图谱。公开矩阵规模不小，细胞也足够多，很容易让人误以为：只要把不同严重程度的细胞放进同一张 UMAP，再跑一轮差异分析，结论就会相当稳。

问题恰恰出在这里。这批数据的统计信息主要来自供者，单个细胞不能充当独立重复。一个供者贡献几千个细胞，并不会自动变成几千个独立样本。健康组、中度组和重度组的供者数都不大，健康组还包含一个来自另一数据集、测序化学也不同的 HC4。只看细胞层面的图和 P 值，很容易把个体差异、样本来源和标签不确定性一并压平。

本文从公开 10x 矩阵重新建立了 QC、降维、聚类、宽谱谱系门控、供者层面的组成和程序分析、raw-count pseudobulk，以及公开 VDJ 的克隆账本。决定结论等级的一步，是把 donor×label 验证失败和单供者稀有 cluster 单独拿出来做敏感性分析。结果并不整齐，却更接近这批数据能承担的证据强度：Severe–Healthy 的 Myeloid 程序和 Myeloid pseudobulk 相对稳定；组成、T/NK 程序和 T pseudobulk 则要明显降级。

论文和数据来源

• 论文：Liao M. et al. Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19. Nature Medicine, 2020。
• DOI：10.1038/s41591-020-0901-9
• PubMed：32398875
• GEO：GSE145926
• 论文指定的额外健康供者 HC4：GSE128033 / GSM3660650。

公开输入是每位供者的 10x 矩阵和疾病供者的 VDJ contig 资源。本文没有从 FASTQ 重新比对，也没有重建作者内部的 barcode correction、ambient RNA 处理或 UMI 生成流程。本文所说的“复现”都应放在这个边界内理解：分析覆盖公开矩阵中的全部合格细胞，但没有逐字节重建作者原始计算环境。

研究设计：先问清楚谁是重复

这项研究的可识别性由供者数决定，而不是由矩阵里有多少 barcode 决定。主链把四位 Healthy、三位 Moderate、六位 Severe 作为 13 个供者单位；其中 HC4 来自 GSE128033，既改变了来源，也改变了 3′ chemistry。HC4 的健康身份和数据来源绑定在同一个样本上，所以不能把“纳入 HC4 后的 Healthy”当作一个没有技术混杂的基线。

这三条链不是三个独立队列。它们共享同一批公开供者，只是改变了外部供者或细胞保留规则。把三条链放在同一张表里，是为了检查一个方向是否依赖特定来源和阈值；不能把重复出现的细胞当作额外的生物学重复。Moderate 只有 3 位供者，组间 exact permutation 的最小可达 P 值也因此较粗，这一点会直接影响组成和程序结果的措辞。

原文的研究问题是描绘 COVID-19 BALF 的免疫细胞景观，并比较健康、中度和重度病例。这里保留这个问题作为背景，却把可检验的问题收紧为：在公开矩阵层面，哪些 Severe–Healthy 方向能经受来源、标签和供者留一敏感性？这样定义问题，是为了让统计单位与可复核输入在分析开始前就固定下来。

先把三个分母分开

主分析纳入 GSE145926 的 12 位供者，并按论文设计加入外部健康供者 HC4，共得到 13 位供者：4 位 Healthy、3 位 Moderate、6 位 Severe。按论文文字阈值重新执行 author-QC 后，主链保留 66,458 个细胞。论文补充表写的是 66,452 个细胞，两者相差 6。这个差异已经在审计表中逐项保留，不能写成“精确重现论文分母”。

另外两条链用于检查来源和 QC 选择是否改变结论：

• GSE145926-only author-QC：63,740 个细胞、12 位供者；
• GSE145926-only strict-QC：53,583 个细胞、12 位供者。

HC4 既来自另一数据集，也使用不同的 3′ chemistry，来源、技术和健康身份完全混在一起。把主链和 GSE-only 链放在一起比较，可以观察结论对 HC4 的依赖，却不能把两条链写成两个相互独立的队列。

图 1 展示的是供者级分母。C142、C144 的公开矩阵线粒体比例偏高，因此严格 QC 会明显改变部分供者的保留量。这里需要同时报告每位供者和每条 pipeline 的细胞数，单独给一个 pooled total 会掩盖样本间差异。

QC 不是一个总开关

逐供者表还记录了输入 barcode、author-QC 细胞、严格 QC 细胞、基因数、UMI 数和线粒体比例。C142 的输入中位线粒体比例约为 88.3%，C144 约为 93.9%；两者在 author-QC 后分别保留 3,411 和 364 个细胞，严格链进一步降到 2,560 和 273。C143 则保留 17,340 个 author-QC 细胞。把这些样本合成一个 pooled QC 百分比，会把样本质量差异隐藏掉。

QC 结果有三种用途：author-QC 固定主分析分母；strict-QC 检查高线粒体和低复杂度细胞是否改变方向；Scrublet 与稀有 cluster 排除检查下游结果对潜在双细胞或供者特异结构的依赖。它们不能被合并成一个“最干净”的单一对象，因为不同阈值回答的是不同问题。每条链都从公开输入重新计算 retention，不从另一条链复制标签或细胞集合。

论文补充表逐供者的总和为 66,452，而按公开矩阵和论文文字阈值重施得到 66,458。差异由逐供者的 -1、+1、+3、+2 等小偏差组成，已经登记在 denominator reconciliation 表中。它足以阻止“分母逐字节一致”的说法，却不足以否定公开矩阵的全量重分析。因此应并列报告两个数字，后续结果采用 66,458。

矩阵层和计算证据

表达对象包含 33,539 个基因特征。用于 UMAP、marker 和 program score 的 X 是归一化后的浮点层；用于 pseudobulk 的 counts layer 保留公开 raw counts。四个 H5AD 对象都经过 CSR payload、obs 总量和基因检出数核对。主链 raw-count 的精确整数和为 462,876,022。若以 float32 对同一批 counts 归约，总和会变为 462,876,032；相差的 10 来自浮点累加误差，因此不能用于整数闭合或分母结论。

主分析完整处理 66,458 个细胞，运行约 2 小时 37 分，峰值内存约 7.8 GiB。这个记录并不能证明生物学结论成立，但能排除“只拿少量细胞或少数基因做展示”的情况。所有定量表都来自完整合格数据；下采样只用于 UMAP 渲染。

H5AD 的 X 与 counts 在这篇文章中承担不同任务。归一化后的 X 适合 PCA、Harmony、UMAP、marker 证据和 program score；raw counts 才能在 donor×label 之后做整数聚合，再计算 pseudobulk。用 X 的浮点总和替代 raw counts，会改变低表达基因的相对权重，也无法留下整数闭合证据。本文按分析环节明确区分所用层；“全量处理”不等于“每个步骤都使用同一层”。

33,539 个特征里，构建 score 的 control pool 明确排除了 120 个预先声明 marker，形成 33,419 个候选 control features；score mapping 有 2,100 行，实际使用 1,847 个不重复 control。这个记录用于防止 marker 被悄悄放回 control pool，也使 score 可以由 source table 独立重算。它仍是同一 assay 内部的构建审计，不是外部表达验证。

HC4 纳入前后，UMAP 能告诉我们什么

图 2 先保留 HC4，给出 all13 结构背景；下面的图 3 沿用相同坐标字段和渲染规则，只把分析对象换成 GSE145926-only，用来检查 HC4 的来源与 chemistry 是否主导可视结构。两张图共用同一分析链，差别只用于来源敏感性审计。

两张 UMAP 用来检查结构背景和来源敏感性。它们可以回答宽谱谱系是否在低维空间中形成可辨认结构，也能直观看到 HC4 纳入后样本分布怎样变化。它们不承担组间定量推断。UMAP 上的距离不能证明轨迹、疾病进展或批次已经被“消除”。

主链得到 22 个 reference clusters，GSE-only 和 strict 链各得到 20 个。独立 A/B 结构重复中，PCA、Harmony、kNN/connectivity、22-cluster partition 和支持标签均一致，ARI 为 1。这个结果说明该计算流程在相同软件、参数和线程设置下可重复，并不等于 22 个 cluster 就是唯一正确的生物分类。

标签证据：cluster 通过，不代表每个 donor×label 都通过

标签构建和验证使用了彼此不重叠的 marker panel。construction markers 负责形成宽谱标签，held-out markers 不进入命名，也不进入 HVG。这样做可以减少“用同一批 marker 命名，再用同一批 marker 证明命名正确”的循环。

图 4 检查标签锁定后的 held-out marker 检出；图 5 随后检查 seed、resolution 与逐供者省略时的标签稳定性。前者回答表达证据是否相容，后者回答计算结构是否稳定，两类证据分别登记，不能互相替代。

在 all13 的 22 个 clusters 中，18 个非 Ambiguous clusters 通过 cluster-level held-out support。GSE-only 为 13/13，strict 为 15/15。问题在于，同一个 broad label 在不同供者里不一定都得到同样的 marker 支持。

把证据下沉到 donor×label 后，一共得到 65 个单元，其中 51 个 PASS、14 个 FAIL，失败单元合计 550 个细胞。组成敏感性使用两种处置：把失败细胞重映射为 Ambiguous，或从相应谱系分母中排除。program 与 pseudobulk 的外部敏感性只重算了谱系排除分支。各分支回答的问题和可用输出不同，不能把其中一条包装成对全部分析层的共同验证。

还存在两个单供者稀有 cluster：cluster 20 是 C152 独有的 33 个 Plasma cells，cluster 21 是 C152 独有的 39 个 Ambiguous cells。按“少于 2 位供者”和“少于 3 位供者”定义，排除集合完全相同，共 72 个细胞。它们可以保留为审计对象，不能被包装成跨供者稳定的新亚型。

从 cluster 到 broad label：证据链的每一层都要留名

标签不是作者原始 taxonomy 的直接复制。本文分析链先用 construction marker 和无监督 cluster 形成 Myeloid、T、NK、Plasma、Epithelial、Ambiguous 六个 broad states，再用不参与构建的 held-out marker 检查 cluster-level support，再把证据下沉到 donor×label。这里的“支持”表示这一 assay 的独立 marker panel 与该标签相容，不表示有外部 atlas 或蛋白实验确认。

标记分数的 control 选择也有审计记录。120 个预声明 marker 被排除在 control pool 之外，HVG 选择没有把 held-out marker 带回建模输入；目标—control mapping、实际使用的 control、匹配距离中位数和 95% 分位数都写入表中。把 marker 画回 UMAP 只能说明构建链内部的一致性，不能把同一张图当成独立验证。

donor×label 的 65 个单元中，51 个通过、14 个失败。失败单元合计 550 个细胞。图 10 同时展示 failure-to-Ambiguous 与 lineage-exclusion 两种组成处置；图 11 和图 12 只来自 lineage-exclusion，因为 failure-to-Ambiguous 分支没有重算 program 或 pseudobulk。这样做把“标签不确定性”变成可计算的分析分支，也把每个分支实际覆盖的结果层写清楚，避免将组成敏感性误读为对 program 和 pseudobulk 的双分支确认。

组成差异：适合描述，不适合写成稳定发现

组成分析回到 donor 层面，使用 CLR effect、exact permutation、bootstrap、within-family BH 和 leave-one-donor-out。all13、GSE-only 和 strict 三条链的 composition family 都没有 BH<0.05；来源与 QC 敏感性主要体现在点估计、区间和方向稳定性，而不是出现一组可直接升级的显著 composition calls。

如果只看图 6 的点估计，某些 NK 或 Plasma 对比似乎有清楚方向。donor×label 失败敏感性改变了这个判断。

为避免编号和证据层交错，图 10 在主链图 7–9 之后集中展示两种 donor×label 失败处置。组成结果在读到该敏感性图之前，只能保留为供者级点估计与不确定性描述。

统计检验怎样回到供者

组成和 program 的 exact permutation 都在供者标签上计算。Severe–Healthy 的 6 对 4 供者组合有 210 个置换，Moderate–Healthy 的 3 对 4 组合有 35 个置换；T/NK 因部分供者没有可估计特征，实际可用组合还会减少。bootstrap 区间描述供者层效应的不确定性，leave-one-donor-out 检查单个供者是否决定方向。family 内的 BH 只在同一组预先登记的测试中执行，不能把不同分析链的 P 值混在一起。

66,458 个细胞是表达和可视化的观察数，6 位 Severe、4 位 Healthy 才是 Severe–Healthy program 的统计单位。当某个 donor 没有足够的某标签时，该 donor 从该 feature 的估计中退出，并在 loo_ineligible_missing_feature_donors 中登记。缺失不是零效应，不能用 0 填补后再扩大样本量。

组成的主链结果也说明了 exact 检验的现实分辨率。Moderate–Healthy 的 NK 点估计为 1.347，exact P 为 0.0571，BH 后为 0.3429；Severe–Healthy 的 NK 点估计为 0.804，exact P 为 0.1238，BH 后为 0.3143。点估计可用于描述方向，有限供者数、multiple testing 和 donor×label 敏感性共同阻止它成为稳定疾病发现。

Myeloid 程序是这批数据里相对稳的一条线

程序分析只在相应谱系的合格 donor 中比较，不把缺少足够细胞的 donor 填成零。Severe–Healthy 的四个 Myeloid programs 在主链、HC4 排除、strict QC、Scrublet sensitivity、construction sensitivity、donor×label exclusion、rare-cluster exclusion 和 donor LOO 中保持方向。

主链中的四个效应为：

这些数值来自 6 位 Severe 和 4 位 Healthy 供者的观察性对比。在这组共享供者的预设分析流程与敏感性范围内，方向保持一致；这不是独立队列复制，也不能改写成因果机制、治疗效应、预后指标或临床决策依据。

图 7 给出主链与 HC4 排除链的 program 结果；图 11 再单独回答 donor×label 失败单元被排除后，Myeloid 与 T/NK 的证据等级怎样变化。

raw-count donor pseudobulk：显著基因也要看 family 和 LOO

Myeloid pseudobulk 从 raw-count layer 按 donor×held-out-supported label 聚合。每个 contrast family 都登记 full test、BH、LOO effect range 和 sign consistency，也保留无法估计的 family，避免结果表只剩下计算成功的部分。

Severe–Healthy 的 BH<0.05 genes 在 all13、GSE-only、strict 三条链中分别为 5,651、4,474、4,258；三条链共同显著且方向一致的基因为 3,572 个。Figure 08 只展示经过共同 gene union 和逐基因 audit 的代表基因，图上正负效应对称呈现。这里的“显著”来自 donor pseudobulk，不来自把数万细胞当作独立样本。

Figure 08 的筛选不是把最漂亮的基因挑出来。figure08_selection_audit.tsv 记录了共同 gene union、测试状态和方向，图中同时保留正向与负向代表。复用这条链时，应先按 donor-level eligibility 建 pseudobulk，再按 family 计算 BH，把跨 pipeline 的 intersection 和 sign consistency 作为稳定性描述；不能先看全细胞差异，再把显著基因汇总到供者层面。

TCR 只能回答疾病供者内部的克隆账本

9 位疾病供者有公开 contig 文件，可以重算 barcode accounting、productive contig、TRA+TRB pairing、clone count 和 Shannon 等描述量。健康供者没有可比 VDJ，因此这里不能形成 Healthy–Moderate–Severe 的推断对比，也不能支持免疫保护、疗效或机制结论。

VDJ 账本还显示每位疾病供者从 contig barcode 到 productive TRA+TRB pairing 的损失链。C141 有 1,065 个输入 barcode，其中 124 个落在 author-QC 且 held-out-supported T 集合之外，852 个 barcode 进入 paired productive 计数；C143 的输入量更大，但同样不能用 barcode 数量替代克隆独立性。multi-chain paired barcode、unpaired productive barcode 和 paired barcode 分开登记，避免把所有 productive contig 当成一个克隆。健康供者缺少同层级 VDJ，任何跨组 TCR 结论都应转化为未来数据需求。

donor×label 失败敏感性：组成、程序和 pseudobulk 分开看

65 个 donor×label 单元中有 14 个 held-out 验证失败。组成层同时计算 failure-to-Ambiguous 与 lineage-exclusion；program 和 pseudobulk 只计算 lineage-exclusion。下面三张图按结果层依次排列，不能把图 10 的双分支设计套到图 11–12 上。

两个组成分支中的 12 个效应都发生数值变化，并且各有 4/12 方向翻转，主要落在 NK 和 Plasma。failure-to-Ambiguous 分支的最大 |effect change| 为 2.0673，lineage-exclusion 分支为 2.0045；前者是图 10 两个分支合并后的最大变化。这个幅度足以改变叙事。组成结果应作为供者级描述和敏感性证据，并明确哪些方向在不同分母下不一致。它不能支撑“重症患者稳定增加或减少某一宽谱谱系”的广泛结论。

图 10 说明组成分母会受失败细胞处置影响。图 11 换到 lineage-restricted program，并且只比较主链与谱系排除分支。

Myeloid Severe–Healthy 的四个 program 方向保持；T/NK 在失败单元排除后，Healthy 的 eligible donors 从 3 降到 2，exact P 的可达分辨率随之变粗。T/NK 结果可以作为后续假设，不能和 Myeloid 放在同一证据等级上。

图 12 继续使用 lineage-exclusion，但统计对象变成 pseudobulk 的 gene×contrast 测试行。同一基因可以在 Moderate–Healthy 和 Severe–Healthy 中各出现一次，因此行数不能写成独立基因数。

Myeloid 有 25,762 个两侧共同测试的 gene×contrast 行，涉及 13,493 个独立基因；没有方向翻转行。active 与 lineage-exclusion 两侧各有 5,651 个 BH<0.05 gene×contrast 行，在该比较表中恰好也对应 5,651 个独立基因。rare-cluster exclusion 同样没有改变 Myeloid 主结果。

T 有 27,742 个共同测试行，涉及 14,688 个独立基因；其中 4,117 个行发生 effect sign flip，涉及 3,511 个独立基因。按 contrast 拆开，Moderate–Healthy 为 13,250 个测试行、1,802 个翻转行，Severe–Healthy 为 14,492 个测试行、2,315 个翻转行。T 在 active 和 lineage-exclusion 两侧都没有 BH<0.05 行。这个结果阻止“所有 tested genes 的方向都稳健”进入结论，也说明行级方向变化不能被误写成 4,117 个互不重复的基因。

三条主 pipeline 与标签敏感性可以合并成一份审计结论：Myeloid 的 BH<0.05 基因数在 all13、GSE-only、strict 中分别为 5,651、4,474、4,258，三链共同显著且方向一致的基因为 3,572；donor×label exclusion 后，Myeloid 的 5,651 个 active call 保持。T 的 27,742 个共同测试行中没有 BH<0.05，且 4,117 个行发生方向翻转。pseudobulk family 没有为每个 LOO 场景重新计算完整的 exact P/BH，这一边界仍需保留。

与原文对照：复现了什么，又缺了什么

本文的公开矩阵分析复核了主要供者分组、论文文字 QC 阈值、宽谱免疫结构，以及重症组 Myeloid 炎症/IFN 上升、alveolar macrophage 和 antigen-presentation 方向下降的信号。与原文相符的地方，需要放在公开矩阵、本文标签门控和供者数的边界内解释。

没有完成的部分同样清楚：没有从 FASTQ 重跑比对，没有作者内部对象，无法恢复 deposition 前的过滤和技术处理；HC4 也不能充当独立健康队列；健康 VDJ 缺失；早期疫情样本的治疗、采样和临床背景不足以支撑当代变异株或临床决策外推。

原文与本文分析的比较可以分成三层。第一层是数据身份和大方向：论文、PMID、GEO accession、供者分组均一致，公开矩阵中 Myeloid 炎症/IFN 程序与 alveolar macrophage、antigen-presentation 方向相容。第二层是可重算但不等价的层：公开矩阵可以重建 QC、宽谱结构、供者 program 和 raw-count pseudobulk，却没有作者内部对象，不能声称相同的 barcode 处理、环境 RNA 修正或对象字段。第三层是公开数据无法承担的层：没有健康 VDJ、没有足够 Moderate 供者、没有同平台 HC4 对照，也没有临床随访，因此不能把关联升级为机制、疗效或预后。

深度解读：这篇复现应该怎么读

这篇分析回答了什么

它回答的是一个较窄的问题：在公开矩阵和现有供者规模下，哪些组间方向经得住 source/QC、标签失败、rare cluster 和 LOO 的连续检查。结果把证据集中到 Severe–Healthy 的 Myeloid program 与 Myeloid donor pseudobulk，其他线索按敏感性强弱分层。

raw counts 和 processed layer 的边界

raw counts 进入 pseudobulk；归一化表达用于结构、marker 和 program。两类层不能混用，也不能用 log-normalized X 的浮点总和替代 raw-count audit。公开矩阵重分析仍然不同于 FASTQ 级重建。

主图应该怎么读

Figure 07 和 08 给出相对稳定的 Myeloid 主链；Figure 10 以两种处置检查 composition，Figure 11–12 只以 lineage-exclusion 检查 program 与 pseudobulk。几张图需要结合各自分支范围阅读。Figure 06 的 composition 点估计如果脱离 Figure 10，会产生过强结论；Figure 09 如果脱离 healthy VDJ 缺失，也会被误读成组间 TCR 结果。

哪些地方不能过度解释

供者数小，HC4 来源与组别混杂，标签仍是宽谱门控，Moderate 只有 3 位供者，T/NK 的 eligible Healthy donors 在敏感性中降到 2。UMAP 不支持轨迹，程序分数不支持功能机制，Myeloid association 也不等于治疗或预后效应。

如果继续深入

更优先的工作，是用独立、同平台或多平台的 COVID/viral pneumonia 队列验证 Myeloid program 和 pseudobulk 方向，同时补充完整临床协变量、健康 VDJ 或正交免疫表型。新的证据才能区分这些方向究竟属于早期疫情 BALF 的队列特征，还是更广泛的严重肺部炎症信号。

对类似旧项目的可复用工作单

这批 BALF 数据适合作为一个完整的审计范例。拿到类似项目时，可以沿着同一顺序建立四张表：输入与分母表把样本、供者、矩阵层和每一步保留数连起来；标签证据表把 construction marker、held-out marker、donor×label 状态和 ambiguous 处置连起来；推断表把生物学重复、测试 family、P 值分辨率、BH 和 LOO 连起来；图表 provenance 表把 source table、脚本、参数、输入对象和最终图片 hash 连起来。

这种组织方式能直接回答旧项目里常见的问题：某个“差异”究竟来自细胞比例还是状态分数？某个 marker 是独立证据还是参与命名的同一批基因？某个显著结果是否由一个供者或一个外部样本驱动？如果这些问题在证据表中没有答案，进一步解释只会掩盖不确定性。

结论审计

当前数据支持

• 13 位供者、66,458 个 author-QC 细胞的公开矩阵重分析，另有 GSE-only 和 strict 两条 sensitivity；
• broad-lineage held-out marker gate、donor×label failure 和 rare-cluster audit；
• Severe–Healthy Myeloid programs 的供者级稳定方向；
• Myeloid raw-count donor pseudobulk 的跨 pipeline、BH 和 LOO 证据；
• 疾病供者内部的 TCR descriptive accounting。

可以作为后续假设

• Moderate、T/NK 和部分 composition 方向；
• Ambiguous/rare states 的进一步参考映射；
• Myeloid program 在独立病毒性肺炎或新冠队列中的验证。

当前不能说

• 与论文完成了相同分母、相同对象或 FASTQ 级完整复现；
• NK、Plasma、T/NK 或全部组成都形成稳定疾病发现；
• Myeloid 信号证明了治疗机制、因果链、预后或临床用途；
• TCR 有健康对照，或可据此判断免疫保护。

这批数据的用处，在于展示供者单位、标签失败和敏感性分析怎样改变结论等级。把这些边界保留下来，公开单细胞矩阵才可能成为可复核的研究证据，避免只留下被细胞数量放大的图。