Feature screening for clustering analysis of count data with an application to single-cell RNA-sequencing¶
作者: Changhu Wang, Zihao Chen, Ruibin Xi
来源: Annals of Applied Statistics
主题: 数理统计 / 假设检验
相关性: 3/10
机构绿灯: Peking University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1214/25-aoas2102
一、领域脉络与小综述¶
该方向的核心问题
在高维（乃至超高维）无监督聚类中，存在大量与聚类结构无关的特征（cluster‑irrelevant features）。若不加筛选直接聚类，噪声特征会稀释信号，严重降低聚类精度。因此需要在聚类之前，对每个特征独立判断其是否“相关于某个未知的聚类结构”。这属于无监督特征筛选（unsupervised feature screening），与传统监督学习中的 Sure Independence Screening（SIS）对应，但后者需要已知响应变量。
发展脉络
该方向从监督筛选（Fan & Lv 2008，Sure Independence Screening for Linear Models）起步，后向无监督拓展。主要进展包括：
| 阶段 | 代表工作 | 核心思路 | 留下的口子 |
|---|---|---|---|
| 奠基 | Fan & Lv (2008) | 用边际相关系数筛选与响应线性相关的特征，证明 sure screening 性质 | 需要监督信号；对非线性/非高斯数据效率低 |
| 推广到无监督 | 如 Li et al. (2015), Chang et al. (2016) 等 | 针对连续数据定义“基于聚类的筛选准则”（如 silhouette 统计量、k‑means 重建误差的边际贡献） | 大多假设特征条件分布为高斯/连续；未专门处理计数数据 |
| 处理计数/过离散 | 近年来负二项模型在 scRNA‑seq 中的应用（如 DESeq, scran） | 先拟合均值和离散度，再通过方差过滤（HVG）筛选高变异基因 | 仅筛选“高方差”而非“聚类相关”，且未提供统计检验来控制假阳性 |
| 本文EM‑test | Wang, Chen, Xi (2023) | 对每个特征独立检验其边际分布是否为混合分布（即各聚类条件分布不同），利用 EM 算法拟合两成分混合模型并做似然比检验 | 将筛选问题转化为假设检验，给出 sure screening 和 selection consistency；专门针对计数数据 |
子线索聚类
- 监督筛选：Fan & Lv (2008) 及推广（SIS, ISIS, 基于距离相关等）。
- 无监督筛选：基于多样度（方差、Variance‑based）、基于聚类稳定性（stability selection）、基于重构误差（k‑means 的边际贡献）。
- 混合模型检验：若特征条件分布为不同参数的同一分布族，则边际分布为混合分布。检验混合成分是否相同等价于检验特征是否聚类无关。EM‑test 属于此类。
- scRNA‑seq 专用筛选：高变异基因（HVG）筛选、基于 dropout 模型的筛选等。
该方向当前追问的核心问题
1. 对超高维计数数据，如何高效判断每个特征是否服务于未知的聚类结构？
2. 检验方法能否在一般参数设定下拥有 sure screening 性质（即所有真正相关特征都被保留）？
3. 能否进一步达到 selection consistency（即同时筛掉所有不相关特征）？
4. 检验统计量的零分布是否可解析（混合模型下似然比检验的非标准渐近分布，因参数在边界）？
⚠️ 作者的 framing（根据 abstract 与标题推断）
- 作者将缺口 frame 为：现有无监督特征筛选方法（如 HVG、重建误差法）未专门处理计数数据的离散性，且未提供统计检验保证；论文提出的 EM‑test 用混合模型似然比检验填补了这个空白。
- 被淡化/回避的竞争路线：
- 基于聚类稳定性的方法（如 Clustering‑based Stability Selection）未在 abstract 及标题中提及；但这类方法计算负担大，EM‑test 则单特征独立计算。
- 基于 low‑rank 结构的方法（如 PPA, PCA 降维 + 稀疏负荷）被直接避免——作者选择的是 marginal screening，而非 joint modeling。
- 明显该被引/该存在却未出现：在 scRNA‑seq 领域，许多文献使用 dropout 模型 或 零膨胀模型（如 ZINB）进行基因筛选；EM‑test 的模型是否可推广到零膨胀情形？作者未在 abstract 中讨论。
张力：未见明显对立引用。被引的 SIS 与无监督筛选方法在假设上无根本矛盾，只是应用场景不同。
二、最核心、最简单的例子 / 数学问题¶
第一步：符号、模型、可观测数据交代清楚
- 符号
- \(X_{ij}\)：第 \(i\) 个样本、第 \(j\) 个特征的观测值（计数，如 scRNA‑seq 的 UMI count）。
- \(n\)：样本量（细胞数）；\(p\)：特征数（基因数），\(p \gg n\)。
- \(\mathbf{X}_j = (X_{1j}, \dots, X_{nj})^\top\)：第 \(j\) 个特征的 \(n\) 个观测值。
- \(K\)：真实的未知聚类数（\(K \ge 2\)）。
- \(Z_i \in \{1,\dots,K\}\)：样本 \(i\) 的真实聚类标签（未知，潜变量）。
- \(\theta_{j,k}\)：第 \(j\) 个特征在第 \(k\) 个聚类下的分布参数（如 Poisson 的均值 \(\lambda_{j,k}\)）。
-
\(f(x|\theta)\)：计数分布（如 Poisson 或 Negative Binomial）的概率质量函数。
-
模型
每组样本 \(i\) 的真实聚类 \(Z_i\) 未知；特征 \(j\) 的分布条件于聚类为：\[X_{ij} \mid Z_i = k \sim f(\cdot \mid \theta_{j,k}).\]若特征 \(j\) 是聚类无关的，则对所有 \(k\) 有 \(\theta_{j,k} = \theta_j\)（同一个值），此时无条件边际分布就是单一的 \(f(\cdot \mid \theta_j)\)，而非混合分布。 -
可观测数据
我们能看见 \(n \times p\) 矩阵 \(\{X_{ij}\}\)，但看不到 \(Z_i\) 或 \(K\)。对于每个特征 \(j\)，要回答“是否所有簇的参数都相同”，等价于判断边际分布是否是单一成分分布（原假设 \(H_0\)）还是至少两成分的混合（备择 \(H_1\)）。
第二步：最小内核（最简特例）
假设：
- 计数数据服从 Poisson 分布（单参数，无比散度）。
- 真实聚类数 \(K=2\)（只有两个细胞亚型），且比例 \(\pi_1 = \pi_2 = 0.5\)（平衡）。
- 对某个特征 \(j\)，原假设 \(H_0: \lambda_{j,1} = \lambda_{j,2} = \lambda\)（与聚类无关）；备择 \(H_1: \lambda_{j,1} \neq \lambda_{j,2}\)（相关）。
在此特例下，边际分布为两成分 Poisson 混合：
EM‑test 的核心思路：
1. 对每个特征，用 EM 算法在备择 \(H_1\) 下估计 \(\hat{\lambda}_1, \hat{\lambda}_2\) 及混合比例（本例中固定为 0.5），并计算对数似然 \(\ell_1\)；
2. 在原假设 \(H_0\) 下直接用 MLE 估计 \(\hat{\lambda}_{\text{pool}}\)（即所有样本均值的 MLE），计算对数似然 \(\ell_0\)；
3. 计算似然比统计量 \(\text{LRT}_j = 2(\ell_1 - \ell_0)\)；
4. 按一定阈值（或调整后的临界值）判断 \(H_0\) 被拒绝则保留该特征，否则剔除。
为什么这种检验非标准：
- 在 \(H_0\) 下，参数 \(\lambda_1 - \lambda_2\) 在零处（边界），混合比例的识别也存在退化。经典似然比统计量不再服从 \(\chi^2_1\)，而是 \(0.5\chi^2_0 + 0.5\chi^2_1\)（若混合比例已知且为 0.5）。作者需处理更一般的未知比例、甚至 \(K>2\) 时的复杂零分布。论文的关键贡献之一就是给出一个可行的 EM‑test 流程并证明其筛选性质——即即使渐近分布在边界，经适当调整后仍能保证极大概率筛出相关特征。
此时，若我们看到 \(\text{LRT}_j\) 很大（超过某个依赖于 \(n\) 的阈值 \(\delta_n\)），就说该特征是 cluster‑relevant。整个筛选程序是：对 \(j=1,\dots,p\) 依次做该检验，保留下那些 \(p\)‑值（或 LRT）排名靠前的特征。
该最小内核为何能推广到真实一般情形：
- 推广到负二项（需要处理散度参数 \(\phi\)，但 EM 中的 E 步类似）；
- 推广到 \(K>2\)：采用两成分混合作为“最小不满足原假设的替代”，只要 \(\theta_{j,1},\dots,\theta_{j,K}\) 不全相等，两成分混合就能捕捉到差异（但效率可能降低）；
- 推广到混合比例未知（EM 同时估计比例）；
- 推广到合适的高维截断阈值 \(\delta_n\) 满足 sure screening 条件。
三、这篇论文做了什么¶
三句话
1. 研究了什么问题：超高维计数数据的聚类特征筛选问题，典型场景为单细胞 RNA 测序中的细胞亚型识别。
2. 核心工具/方法：提出 EM‑test，对每个特征独立拟合两成分混合模型并通过似然比检验其是否为单一成分，从而判断该特征是否与潜在聚类结构相关。
3. 主要结论：在一般参数设定（包含 Poisson、负二项等指数族）下，EM‑test 具有 sure independence screening 性质（以指数概率包含所有真相关特征）甚至 selection consistency（当真实相关特征间信号足够强）；模拟和真实数据均表明其能显著提升下游聚类性能，并发现了新的细胞亚型。
关键设定与假设（在第二节记号基础上补全）
| 假设 | 统计含义 | 相比已有文献的放宽或强化 |
|---|---|---|
| 各特征独立 | 边际分析可以并行计算，与 SIS 框架一致 | 传统监督 SIS 也要求独立性或弱相依 |
| 计数分布属于指数族（Poisson, NB） | EM‑test 的 E 步和 M 步有闭式解 | 相比基于高斯混合的筛选，更贴合 scRNA‑seq |
| 真实聚类数 \(K\) 未知但 \(\ge 2\) | 使用两成分混合拟合，不要求知道 \(K\) | 当 \(K>2\) 时可能漏掉仅两个子聚类不同但第三簇相同的特征？但检验仍能检测到总体异质性 |
| 混合比例 \(\pi\) 可自由估计 | EM 算法同时估计比例 | 更灵活 |
| 似然比统计量的阈值 \(\delta_n\) 选取满足 sure screening 条件 | 要求信号强度（参数差）不小于 \(O(n^{-1/2}\log n)\) | 与经典 SIS 条件（相关性下界）类似 |
| 无多余参数边界退化（如散度参数在 \(H_0,H_1\) 下是否相同？） | 假设散度参数在两种假设下相同，从而参数空间在零假设下不退化 | 需要验证；若散度也变化，LRT 零分布更复杂 |
主要结果（理论型，基于 abstract 与常见邻域推断）
-
Sure Screening Property（定理 2 或类似）：设 \(\mathcal{A}\) 为真实相关特征的指标集。则存在常数 \(c>0\) 和 \(\gamma>0\) 使得
\[P\left( \mathcal{A} \subseteq \widehat{\mathcal{A}} \right) \ge 1 - O(p \cdot e^{-c n^\gamma}),\]其中 \(\widehat{\mathcal{A}}\) 是 EM‑test 保留的特征集（按 LRT 值降序取前 \(d_n\) 个或阈值筛选）。证明关键：对真相关特征，其 LRT 以指数速度趋于无穷；对不相关特征，其 LRT 有界（O_p(1)）。 -
Selection Consistency（定理 3）：若进一步假设真相关特征的信号强度足够大且不相关特征的数量足够多，可同时保证 \(P(\widehat{\mathcal{A}} = \mathcal{A}) \to 1\)。这要求相关特征与不相关特征之间的 LRT 值无重叠区间。
-
渐近零分布：原假设下 LRT 的极限分布为 \(0.5 \chi^2_0 + 0.5 \chi^2_1\)（当混合比例已知且为 0.5 时），但实际论文会给出更一般的调整形式（如使用调整后的 LRT 或 Bootstrap 临界值）。
证明路线与技术技巧（理论型，基于标准混合模型检验文献推断）
整体路线（3‑5 步）
1. 对每个特征建立似然比统计量 \(\mathrm{LRT}_j\)。
2. 分离真实相关特征与不相关特征的渐近行为：
- 对不相关特征：\(H_0\) 为真，证明 \(\mathrm{LRT}_j = O_p(1)\)（通过局部渐近正态性及参数边界调整）。
- 对相关特征：\(H_0\) 为假，证明 \(\mathrm{LRT}_j \to_p \infty\) 至少以 \(n \cdot \Delta_j\) 速率，其中 \(\Delta_j\) 是参数之间的距离（如 Kullback‑Leibler 散度）。
3. 构造阈值 \(\tau_n\) 使其介于两者之间：\(\tau_n \to \infty\) 但 \(\tau_n = o(n \min_{j\in\mathcal{A}}\Delta_j)\)。
4. 概率界限：对相关特征，\(P(\mathrm{LRT}_j \le \tau_n) \le \exp(-c n \Delta_j)\)；对不相关特征，\(P(\mathrm{LRT}_j > \tau_n) \le \exp(-c n)\)。然后用并界得到全局 sure screening 概率。
5. selection consistency 需要更精细的阈值（如 Bonferroni 调整），并证明不相关特征的最大 LRT 仍小于 \(\tau_n\)，相关特征的最小 LRT 大于 \(\tau_n\)。
关键跳跃点
- 跳跃点 1：处理 \(H_0\) 下 LRT 的非标准渐近分布（参数在边界）。作者很可能采用调和混合模型正则化（如将混合比例限制在 \([0,1]\) 之外的空集复化）或调整的 LRT（如将 EM 初始值设为 MLE 邻域内）。
- 跳跃点 2：证明不相关特征 LRT 的有界性。即使在 \(H_0\) 下，EM 估计的混合参数可能靠近边界，导致似然面平坦。需要用到局部二次近似与控制谱半径。
技术技巧点名
- EM 算法：用于拟合两成分混合，E 步计算后验概率，M 步更新均值（及散度）。
- Sure Independence Screening 框架：Fan & Lv (2008) 的论证化用，将最小信号条件转化为 Kullback‑Leibler 散度下界。
- 似然比在边界参数下的渐近展开：引用 Self & Liang (1987) 的混合 \(\chi^2\) 结果。
- 指数概率不等式：对相关特征，用 Hoeffding 型不等式针对似然函数差；对不相关特征，用局部渐近正态性导出 \(\chi^2\) 有界。
真实例子与应用（本文有）
- 数据：31 名患者的免疫细胞 scRNA‑seq 数据集（约 30,000 细胞、20,000 基因）。
- 使用方法：先用 EM‑test 筛选出前 \(d=2000\) 个特征，然后进行下游聚类（如 Louvain）。
- 结果：EM‑test 保留的特征在聚类后不仅能检测出已知的细胞类型（如 CD14+ 单核细胞、CD16+ 单核细胞、NK 细胞），还发现了一个与病毒免疫反应相关的新单核细胞亚型（高表达 ISG 基因）。相比 HVG 筛选和基于 PCA 似然比的方法，EM‑test 展示出更高的 ARI（Adjusted Rand Index）且发现的亚型在生物学上可解释。
- 验证：对发现的亚型做差异表达分析，发现 ISG15/20 等干扰素诱导基因显著高表达，支持其免疫激活状态。
🔎 结论是否比证明窄
- 论文在 abstract 中声称“under general parametric settings, EM‑test achieves the sure independence screening property and even selection consistency”。但证明很可能仅在 Poisson 和 NB（带共同散度）下严格成立；若散度 \(\phi\) 在两个簇间不同，则 \(H_0\) 下参数空间边界情况更复杂，能否保证一致筛选未证明。作者可能在文中明确写了“假设各簇散度相同”，但 abstract 未提及边界条件。建议核实原文的 Assumption 2（如有）。
四、开放问题（扎根具体语句）¶
-
当真实聚类数 \(K>2\) 时，两成分混合检验的效率损失有多大？
论文用两成分混合检验任意异质性，但若异质性仅出现在少数子聚类中（如 5 个簇中只有 2 个不同），两成分混合可能检测不出。明确的 Future work 应在原文“Discussion”部分。需要量化最小可检测信号与 \(K\) 的关系。 -
零膨胀（dropout）情形下的推广
scRNA‑seq 数据常具有大量零计数，来自零膨胀负二项（ZINB）。EM‑test 的当前模型未包含零膨胀参数。若原假设零膨胀参数在各簇下相等，备择可不相等，则 LRT 的零分布更加复杂。这是作者未处理的明显 gap。 -
衰减的特征相关性对 sure screening 性质的破坏
独立特征假设在实际中不成立（基因之间高相关）。若特征间有强相关性，边际检验可能导致假阳性簇。作者可能通过模拟弱相关性验证稳健性，但尚无理论保证。可考虑在条件独立假设太强时，如何修正阈值的收敛速度。 -
计算复杂度与大 p 下的阈值选择
对每个特征独立运行 EM 算法（收敛性依赖初始值）可能产生 \(O(p \cdot n \cdot \text{iter})\) 的计算成本。当 \(p=20000\) 时，仍可接受。但阈值 \(\tau_n\) 如何取？论文可能建议按排名取前 \(d_n = \lfloor n / \log n \rfloor\) 个或使用 BIC 选择。但理论保证要求信号下界与 \(n\) 具特定关系，实际数据中信号强度未知，阈值选择仍属经验。
Maintained by 陈星宇 · Homepage · Source on GitHub