论文解析-CCST-基于GNN的空间转录组聚类算法

参考

Li, J., Chen, S., Pan, X. et al. Cell clustering for spatial transcriptomics data with graph neural networks. Nat Comput Sci 2, 399–408 (2022). https://doi.org/10.1038/s43588-022-00266-5

动机

1. 在体外培养的相同细胞周期的细胞更有可能在一起
2. 体内组织的某些细胞类型已知在空间上与自身或特定的细胞类型距离近
   因此，可以利用空间结构作为改进细胞聚类的有效信息。已有方法大多假设相同组织的细胞空间上相近，但没有考虑组织样本中复杂的全局细胞相互作用。

数据集

Dataset	Genes	Cells	Cell Types	Others
MERFISH	10,050	1368	-	human osteosarcoma cells from 3 batches
seqFISH+	10,000	2.050	11	mouse OB (cortex)
DLPFC	-	-	12 samples,up to six cortical layers and white matter	human dorsolateral prefrontal cortex
10× Visium spatial transcriptomics	-	-	human breast cancer

方法

DGI

已有的一种以无监督方式学习图结构数据中节点表示的方法，DGI依赖于最大化训练后的节点表示与图整体特征的互信息

本文通过Readout函数得到整图的特征，即通过聚合节点特征的方式来得到整图的特征表示。最常用的Readout函数为
$$R\left(H\right)=\sigma \left(\frac{1}{N}\sum _{i=1}^N{\overrightarrow{h}}_i\right)$$
即整图特征=激活函数（每个节点平均特征）
更多DGI相关信息可参考这里
更多Readout相关信息可参考这里

混合图邻接矩阵

1. 构建单一批次（batch）的邻接矩阵
   计算细胞i和j的距离dij，如果大于阈值则认为有边，如果小于阈值则认为无边
2. 跨批次细胞
   如果细胞i和j不在同一batch中，则认为无边
3. 混合图邻接矩阵
   
   其中I表示单位矩阵。λ为超参数，空间转录组数据中为0.3，单细胞转录组数据中为0.8

数据预处理

单细胞转录组数据

1. 移除低质量基因
   删除在所有细胞中不表达的基因
2. 校正批效应
   Scanorama算法
3. 表达水平标准化
   
   其中，i表示细胞，j表示基因
4. 删除低变异基因
   用标准化后的基因表达水平计算基因方差，删除方差低于0.4的基因。

空间转录组数据

1. 删除低质基因
   （1）删除数据集中表达少于3个位点的基因
   （2）删除每个细胞中平均表达水平低于0.02的基因
   （3）删除方差低于0.05的基因
2. 表达水平标准化
   同上

图结构

DGI

无监督图嵌入方法DGI，形成从原始输入到隐变量的映射。即学习一个映射函数 E ( X , A ) = H = { h 1 ,   h 2 ,   . . . ,   h N } E\left( X,A \right) =H=\left\{ h_1,\ h_2,\ ...,\ h_N \right\} E(X,A)=H={h1​, h2​, ..., hN​}
其中X为原始节点特征，A为邻接矩阵，H为训练后的节点特征

4层GCN

（1）GCN公式
本文采用标准GCN

（2）激活函数
本文采用参数整流线性单位函数作为激活函数，a为可学习的参数

（3）生成负样本邻接矩阵
$$\bar{A}=C\left(A\right)$$
即负样本邻接矩阵$\bar{A}$为随机图结构，因此可以根据该随机拓扑结构得到每个节点的局部表征${\bar{h}}_i$
（4）目标函数
目标函数为DGI目标函数，最大化L，即最大化正样本节点与整图特征的互信息，最小化负样本节点与整图特征的互信息

D为互信息，E表示由不同网络结构训练的模型，s为由正样本图结构训练的得到的整图特征（利用Readout函数由每个节点得到整图特征）

聚类

1. DGI输出每个节点训练后的特征，用PCA进行降维
2. 在第一主成分维度上用k-means++（在选取初始种子节点时，离已选种子节点越远的节点被选中的概率越大参考）进行聚类
3. 如果未知类别数，则根据轮廓系数（silhouette score）确定聚类数

结果

应用CCST于空间基因表达数据

数据集：带有3个复制本的MERFISH，每一列对应一个复制本
第一行为细胞的空间位置，不同颜色对应不同聚类
第二行为该复制本构建的图结构，可看出空间位置近的细胞网络结构紧密

计算邻居富集率：对于一个类别的所有细胞，根据图结构收集每个细胞的邻居细胞，计算邻居细胞属于每个类别的占比，得到邻居富集率
结果表明每个类别的邻居富集率均很高
与MERFISH数据集文献对比，每个聚类仅对应一个细胞周期阶段，且第一列中C0与C3空间位置相近，所以邻居富集度较高。

识别细胞周期阶段

1. 在每个聚类中，用Mann–Whitney U test识别高表达差异基因
2. 在整个数据集中选取top200差异表达基因，进行GO富集分析
3. 富集分析结果根据 -log(P value) 升序排序
   
   结果表明，每个聚类对应一个细胞周期阶段，框中为重要GO
   聚类C1中差异表达基因主要对应基因复制功能，说明C1对应S阶段；
   聚类C3中差异表达基因主要对应RNA处理功能，说明C3对应G2阶段
   聚类C0中差异表达基因主要对应有丝分裂功能，说明C0对应M阶段；
   聚类C4中差异表达基因主要对应许多过程的负调控功能，说明C4对应G1阶段
   
   作者定义每个细胞阶段的关键GO，计算每个算法能够识别的关键GO数量。从e图看出CCST识别出的关键GO数量最多，f图说明CCST识别的统计学显著性最高。g图为每个算法对4个聚类识别的GO重复度，CCST重复度低，说明对生物功能分析的显著性和可靠性。

对比SSCT与其他处理ST数据算法的性能

作者对比CCST与7种其他算法在DLPFC数据集上的聚类性能，常用聚类指标ARI、FMI、NMI值越大表明算法越准确。
空间位置相近的细胞通常属于同一大脑皮层，但有些相近的细胞属于不同大脑皮层。为了测试算法对该现象的识别性，本文采用LISI指标（the local inverse Simpson’s index）衡量细胞的局部多样性，值越小表明不同聚类的空间分离性越好。

在seqFISH+数据集上识别细胞亚型

在已被良好标注的数据集上，对每种细胞类型进行二分类，以发现潜在细胞亚型。
a图为CCST对中间神经元细胞发现的细胞亚型
b图为2种细胞亚型的邻居富集度分布（图中隐去中间神经元细胞），说明这2种细胞亚型的微环境差异明显
选取top200差异表达基因进行GO分析，c图对应聚类C0，发现主要是神经功能相关；d图对应C1，发现主要是定位调控和细胞迁移等功能。
因此，该人工注释的中间神经元细胞可以分为2个亚型，一个为可以与其他神经元交流的成熟功能性神经细胞，一个为与邻居细胞交互的未成熟细胞，主要进行定位和迁移。这个结论也可从其他文献中得到验证。
同时，作者还在其他3种人工注释的细胞中识别出新的细胞亚型。