基于常数分段变点数据聚类的贝叶斯方法

《Stats》：A Bayesian Approach for Clustering Constant-Wise Change-Point Data Ana Carolina da Cruz and Camila P. E. de Souza

　　

变点模型用于处理有序数据序列，其核心目标是推断数据序列中某一方面发生变化的位点。研究人员提出并实现了一种非参数贝叶斯模型，用于基于常数分段变点轮廓对观测值进行聚类，通过吉布斯采样器（Gibbs sampler）实现。该模型在常数分段变点结构上引入狄利克雷过程（Dirichlet process, DP），以在聚类的同时进行多重变点估计。此外，该方法能够自动模型中的聚类数量，无需事先指定聚类数目。研究人员在多种模拟场景和单细胞基因组测序的真实数据集上对该方法进行了评估，均获得了令人满意的聚类和估计结果。所提出的方法已实现为名为BayesCPclust的R包，可从综合R档案网络（CRAN）获取。

## 研究背景与问题

变点模型致力于分析有序随机量序列，其应用涵盖从日常平均气温到基因组测序数据等多种场景。变点检测旨在推断数据序列中如位置或分布等某一方面发生变化的位置，具有重要实际意义。变点模型最早可追溯至20世纪50年代，Page等人提出的累积和控制图（CUSUM）用于检测质量控制过程中均值的变化。此后，变点检测在金融分析、生物统计学和时间序列分析等领域得到广泛研究。变点模型通常分为在线方法（online methods）和离线方法（offline methods），前者用于新数据持续到达时的序贯检测，后者则对完整观测序列进行回顾性分析。此外，变点模型可分为参数模型和非参数模型：参数模型假设潜在分布属于某个已知分布族，而非参数方法则依赖于密度函数的估计，但适用范围更广。

现有变点检测方法多针对单一数据序列，然而在多序列间识别共同模式以将来自同一分布的序列进行分组具有重要意义。尽管已有研究尝试结合变点估计和基于模型的聚类，但关于多序列变点数据聚类的研究仍显不足。尤其在单细胞拷贝数数据分析中，现有方法通常采用隐马尔可夫模型（HMMs）推断拷贝数状态后再进行聚类，缺乏能够同时进行拷贝数轮廓分析和聚类的替代方法。

## 研究目的与方法

为此，研究人员开发了一种非参数贝叶斯模型，通过吉布斯采样器对多重常数分段变点数据进行聚类。该模型在常数分段变点结构上引入功能性狄利克雷过程，自动模型中的聚类数量，无需预先指定聚类数目。这是首个为多重变点数据聚类与同时变点检测问题提供实现的工作。该方法应用于基于拷贝数轮廓对异常（肿瘤）单细胞基因组数据进行聚类，并在多种模拟场景下评估了性能。研究成果以R包BayesCPclust的形式实现并发布于CRAN。

该研究的核心技术方法包括：（1）**非参数贝叶克雷过程先验**：在常数分段变点结构上设定狄利克雷过程先验，使模型能够自动确定聚类数量；（2）**吉布斯采样器**：构建贝叶斯推断框架，通过吉布斯采样迭代更新变点结构、截距参数、方差分量、聚类分配及相关超参数；（3）**分层基线分布设计**：基线分布G₀分层定义为：变点数量K服从截断泊松分布，变点间区间长度服从多项分布，常数水平α_l服从?上的非信息先验；（4）**Metropolis-Hastings算法**：用于超参数λ的更新；（5）**后验推断与标签处理**：通过后验众数确定离散参数，采用最频繁标签分配解决标签切换问题。样本队列为NordenSttrom等人发表的单细胞拷贝数数据集（SRP074289）中的结直肠癌患者CRC2的45个细胞。

## 研究结果

**3.1 吉布斯采样器初始化与实现**

研究人员采用两条不同初始值的链进行收敛诊断。第一链从真实设定初始化，第二链则在真实值基础上施加微小扰动。模拟在Sharcnet的Graham集群上完成，使用R 4.2.1版本及parallel、extraDistr、RcppAlgos、MASS和FDRSeg等包实现并行计算。

**3.2 性能评估指标**

每条链运行5000次迭代，舍弃50%作为burn-in期，并对剩余样本进行25倍稀释，最终保留200个样本进行后验推断。截距估计采用后验均值，离散变量采用后验众数。聚类性能采用V度量（V measure）评估，该指标同时衡量聚类的同质性和完整性。

**3.3 情景一：变化数据序列数量（σ_n²≈0.05）**

在该情景中，研究人员在固定其他参数的条件下，变化数据序列数量N=10,25,50。结果显示，该方法能够准确恢复各聚类的真实截距值，且随着数据序列数量增加，截距估计的后验方差和95%可信区间平均长度均减小，表明不确定性降低。所有96个数据集均正确估计了变点位置和聚类分配，V度量均为1。方差分量的平均绝对偏差（MAD）保持稳定，表明增加数据序列数量对提高方差分量估计精度作用有限。

**3.4 情景二：变化数据序列数量（σ_n²≈0.50）**

该情景评估了更高数据离散程度下的方法性能，方差分量均值较情景一增大10倍。结果显示，尽管后验方差和可信区间长度因数据离散度增加而明显增大，但截距估计仍然准确，且随着数据序列增加不确定性减小。变点位置和聚类配置在所有数据集中均被正确恢复，V度量均为1。方差分量估计的MAD较小且稳定，进一步验证了增加数据序列数量对方差分量估计精度提升有限。

**3.5 情景三：变化位置数量**

研究人员变化位置数量M=50,100,200，固定其他参数。结果显示，该方法正确估计了各聚类的截距，且随着位置数量增加，截距估计的不确定性降低。所有数据集均正确识别了变点数量和位置，V度量均为1。值得注意的是，增加位置数量显著提高了方差分量的估计精度，MAD随M增加而降低，这与数据序列数量的影响形成对比，因为每个数据序列具有独立的方差分量。

**4. 真实数据分析**

研究人员将方法应用于结直肠癌患者CRC2的单细胞拷贝数数据，该数据包含45个来自原发肿瘤和肝转移灶的细胞，覆盖染色体19、20、21的583个基因组窗口。为降低计算成本，采用窗口大小为5的中位数滑动窗口将位置缩减至290个，并设定最大变点数量k^*为2。

该方法识别出三个聚类：聚类1包含18个原发肿瘤细胞，在变点位置100和226处有明显变化；聚类2包含4个原发肿瘤和5个转移肿瘤细胞，所有位置的log2比值约等于1，截距不显著；聚类3包含18个转移肿瘤细胞，变点位于165和215。这一结果与NordenSttrom等人的研究部分一致，原发和转移肿瘤细胞基本分开聚类，但亦发现了一个新的混合来源细胞聚类，为数据集提供了新见解。

## 讨论与结论

**6. 结论**

模拟研究表明，该方法能够恢复数据序列的真实分类，且在变化数据序列数量和位置数量时均能精确识别变点。重要的是，数据离散程度不影响方法性能——在方差分量从小平均值的逆伽马分布和大平均值的逆伽马分布中采样时均获得满意结果。此外，该方法对异常值具有鲁棒性：在聚类1的9%序列中将第19个观测值减少10个单位导致变点偏移一个位置时，方法仍能成功恢复真实变点轮廓和正确聚类分配。

在单细胞拷贝数数据集上的应用结果显示，该方法与NordenSttrom等人的观察一致，基于变点结构获得了相似的肿瘤细胞聚类，部分细胞按组织来源聚类。同时，研究发现了一个由原发和转移肿瘤细胞组成的新聚类，为理解肿瘤异质性提供了新视角。

为进一步促进方法的应用，研究人员开发了R包BayesCPclust并发布于CRAN，这是首个同时解决多重变点数据聚类和变点检测问题的软件包。

该方法的局限性在于计算成本：需要计算变点间区间长度所有可能组合的概率，当位置数量增加时计算量庞大。在真实数据分析中，研究人员通过抽样部分可能组合来缓解这一问题，但效果有限。随着数据序列或位置数量增加，处理时间显著延长，最高位置数量场景平均需20-30小时。此外，与MCMC方法一样，标签切换是常见问题，研究人员通过采用最频繁标签集进行校正后获得聚类特定的参数后验估计。