(PNAS 2019) scHiCluster（Part II：

2021-07-26 本文已影响0人阿狸的窝

上一篇：(PNAS 2019) scHiCluster （Part I：文章梳理）

原文：
Zhou J, Ma J, Chen Y, Cheng C, Bao B, Peng J, Sejnowski TJ, Dixon JR, Ecker JR. Robust single-cell Hi-C clustering by convolution- and random-walk-based imputation. Proc Natl Acad Sci U S A. 2019 Jul 9;116(28):14011-14018. doi: 10.1073/pnas.1901423116. Epub 2019 Jun 24. PMID: 31235599; PMCID: PMC6628819.

基于 bulk Hi-C 数据生成单细胞Hi-C模拟数据

作者在文章中指出，单纯地将bulk Hi-C数据进行简单抽样（downsample至相同的total contacts number）不能很好地模拟single-cell Hi-C数据。通过简单抽样得到的结果稀疏度（sparsity）和变异度（variabielity）均更低。

Figure S2 A-C
为解决此问题，作者首先将抽样控制的对象从total contacts变为sparsity。
基于 Ramani et al.[3] 和 Flyamer et al.[4] 两个数据集，作者观察到 total contact $C$ 与sparsity $S$ 间存在对数线性关系：
$log S = a log C + b$

具体步骤

原始数据获取：从Juicebox中提取100 kb resolution下的MAPQ30 contact matrices
分辨率整合：整合得到200 kb 和 1Mb resolution的contact matrices（记作 $A$ ）
归一化校正：使用 SQRTVC normalization

SQRTVC normalization
$D_{ii}= \sum_{j}A_{ij} \\ B = D^{-{1\over 2}} A D^{-{1\over 2}} \\$
稀疏度控制：对于每个单细胞随机产生 $t \sim Union(logM-0.5, logM+0.5)$ ，则该细胞的total contacts为 $C =e^t$ ，根据 $log S = a log C + b$ 计算得到 $S$ ；
添加噪音：当contact matrix包含 $n$ 个bins时，随机生成一个长度为 $n$ 的向量 $R$ ，其中 $R_i \sim Union(-k,k)$ ， $k$ 为 noise level。
$E_{ij} = B_{ij} \times R_{|j-i|}$

6.生成模拟数据：从 $E$ 中选取 $S$ 个位点，将 $C$ 个contacts分配在这些位置

从 $E$ 中随机选取 $S$ 个位点，

作者基于Stevens et al.数据（8 G1-phase mESC），将8个单细胞数据合并到一起产生 pseudobulk 数据集，然后再从pseudobulk数据集中使用4种不同的方式抽样。将抽样产生的模拟scHiC数据与原有的真实scHiC数据放在一起进行主成分分析，作者证明了contol sparsity + noise 组可以更好的模拟真实的scHiC图谱。

截屏2021-07-26 下午5.51.59.png

scHiCluster

Convolution-based imputation

作者认为，scHiC的contact matrix中的0值并不意味着两个位点间不存在interaction，而应当视作由实验的技术手段导致的缺失值，因此需要进行imputation。
作者首先假设，如果位点A和位点B在空间上距离很近，则A应当与位点B在线性基因组上的邻居也很近。因此，作者使用Convolution-based imputation使信息在一维邻居间传递。

给定 window size $w$ 和一个 $m \times m$ 的filter $F$
对于一个 $n \times n$ 的contact matrix $A$ ，imputed matrix
$B_{ij} = \sum_{i-w \leq p \leq i+w\\j-w \leq q \leq j+w} F_{pq}A_{pq}$

默认参数：

$F_{ij}=1$
对于 1 Mb resolution maps， $w=1$

Random-walk-based imputation

Random walk with restarts（RWR） [6-7]
假设一个包含 $N$ 个节点的图 $G$ 具有邻接矩阵 $W$
初始状态下，一名漫步者位于某一初始节点 $v_i$ 。
在每一轮的移动中，walker可以选择移动到任意的与当前节点相连的节点，到达每个点的概率与边权成正比。另外在每一轮漫步者有 $\gamma \in (0,1)$ 的概率返回起点。
$P^{(t)}$ 是一个 $N\times N$ 的矩阵，其中 $P_{ij}^{(t)}$ 为漫步者从 $v_i$ 出发，经过 $t$ 轮移动后位于 $v_j$ 的概率
根据移动规则
$P^{(t+1)}=(1-\gamma) W'_{ij}P^{(t)} + \gamma P^{(0)}$
其中
$W'_{ij}=\frac{W_{ij}}{\sum_jW_{ij}}$
随着 $t$ 的增大， $P(t)$ 会趋于收敛

使用convolution-based imputation步骤得到的 $B$ 作为邻接矩阵，则转移概率矩阵
$C_{ij} = \frac{ B_{ij} }{ \sum_{j}B_{ij} }$
初始化 $Q^{(0)}=I$ (单位阵）
设定restart概率为 $p$ ，迭代计算 $Q^{(t)} = (1-p)Q^{(t-1)}C + pI$
当 $|| Q^{(t)} - Q^{(t-1)} ||_2 < 10^{-6}$ 时，视作收敛，迭代结束

Select top interactions

记RWR收敛于矩阵 $Q$
使用 $Q$ 的80th quantile作为阈值 $t$
将 $Q$ 转化为0/1矩阵 $Q_b$ （即仅选取top 20%的interaction）

Embedding

将 $n \times n$ 的矩阵 $Q_b$ 转化为长度为 $n^2$ 的向量
所有 $m$ 个细胞合并在一起得到 $m \times n^2$ 的矩阵
首先对于每条染色体分别使用PCA进行降维
将各染色体的结果合并到一起，使用PCA进行二次降维
Figure S10中所使用的weight matrix 为两次PCA的weight matrix的点积（dot product)

clustering

使用前10个主成分（PCs）进行K-means ++ 无监督聚类（根据已知的细胞种类数目设定K）

K-means

初始时，随机选取k个样本作为聚类中心；

每一轮计算每个样本到各聚类中心的距离，并将样本分配到预期距离最近的中心所在类

对于每一类，计算该类中所有样本的质心作为该类新的中心

重复步骤2-3直至满足终止条件（小于最小误差或达到最大迭代步数）

K-means++
K-means++在Kmeans的基础上，对初始聚类中心的选择方法进行优化，聚类中心选取方法为：
初始时算法随机选取一个中心点 $a_1$
假设已经选取了 $n$ 个聚类中心，对于剩余的样本，计算每个样本距离已有的 $n$ 个中心的距离，假设其中最远的距离为 $D(x)$ ，则该样本被选为第 $n+1$ 个聚类中心的概率为
$\frac{ D(x)^2 }{ \sum{D(x)^2} }$
即与已有中心距离越远的样本，越有可能被选择为新的聚类中心

聚类效果评价

当细胞真实的所属类别已知时，使用ARI评价聚类的准确性。

Rand Index
假设：

样本 $i$ 的真实类别为 $X_i$ ，聚类时被判定为第 $Y_i$ 类；

样本 $j$ 的真实类别为 $X_j$ ，聚类时被判定为地 $Y_j$ 类

定义：

TP : $X_i =X_j$ 且 $Y_i = Y_j$ （两个样本来自同一类，在聚类结果中也被归入同一类）

FP： $X_i != X_j$ 但 $Y_i = Y_j$ （两个样本属于不同类，但是在聚类结果中被归入同一类）

FN： $X_i = X_j$ 但 $Y_i != Y_j$ （两个样本属于同一类，但是在聚类中被归入不同类）

TN： $X_i != X_j$ 且 $Y_i != Y_j$ （两个样本来自不同类，在聚类中也被归入不同类）

定义
$Rand\ Index = \frac{TP+TN}{TP+FP+TN+FN}$

Adjusted Rand Index
对于两个随机划分，Rand Index的期望不为0，为解决这一问题，定义校正兰德系数
$ARI = \frac{ RI - E(RI) }{ max(RI) - E(RI) }$
ARI的取值在［－1，1］之间，负数代表结果不好，越接近于1越好对任意数量的聚类中心和样本数，随机聚类的ARI都非常接近于0

其他聚类方法

PCA

对于每个细胞，对raw contact matrices进行log2转化
将 $n\times n$ 的矩阵转化为长度为 $n^2$ 的向量
使用与scHiCluster相同的方法进行两轮PCA降维

HiCRep + MDS

使用 1-Mb resolution ，对raw contact matrices进行log2转化
smoothed with window size 1
计算任2个细胞间对于每条染色体的 stratum-adjusted correlation (SCC) 系数，所有染色体的SCC的中位数作为最终的SCC distance（ $d_{scc}$ ）
将SCC距离转化为欧氏距离： $d_{euc}=\sqrt{2-2d_{d_{scc}}}$
基于Euclidean距离矩阵进行MDS降维

Eigenvector

1.对于每个细胞，对raw contact matrices进行log2转化

计算 distance-normalized matrix B
$B_{ij}= \frac{ A_{ij} }{ \sum_{i'} A_{i' + j-i} }$
基于 distance-normalized matrix $B$ 进行PCA降维，提取PC1作为细胞特征
（特殊处理：分别计算 PC1>0 和 PC1<0 的bins的的平均CpG，如果后者更高则取PC1的相反数作为细胞特征）
联合m个细胞得到 $m \times n$ 的特征矩阵，使用PCA进行降维

Decay

此方法提取的是 P(d) ~ d 曲线特征，或者说是距离为 $d$ 的contact占比

对于每个细胞，对raw contact matrices进行log2转化
计算decay 特征
$B_d = \frac{ \sum_j A_{j,j+d} }{ \sum_{i,j}A_{i,j} }$

TAD-like structure identification

TAD in bulk HiC (ESCs and NPCs)：使用TopDom（10 kb resolution, window size=5）计算得到TAD
Figure 5D 中包含了 nondiagonal contacts > 100k 的 1007 个 mESC细胞，使用TopDom（40 kb resolution, window size=5）计算得到TLS

比较 TAD与TLS边界

假设在bulk Hi-C中有TAD边界 $(i,j)$ ，在scHiC中有TLS边界 $(i',j')$
当满足一下条件时，认为 TAD与TLS重合

$|i'-i| \leq 80 kb$ 且 $|i'-i| \leq 0.25 \times (j-i)$
$|j-j'| \leq 80 kb$ 且 $|j'-j | \leq 0.25 \times (j-i)$

参考文献

[1] S. S. P. Rao et al., A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 1665–1680 (2014).
[2] B. Bonev et al., Multiscale 3D genome rewiring during mouse neural development. Cell 171, 557–572.e24 (2017)
[3] T. J. Stevens et al., 3D structures of individual mammalian genomes studied by single- cell Hi-C. Nature 544,59–64 (2017).
[4] V. Ramani et al., Massively multiplex single-cell Hi-C. Nat. Methods 14,263–266 (2017)
[5] I. M. Flyamer et al., Single-nucleus Hi-C reveals unique chromatin reorganization at oocyte-to-zygote transition. Nature 544, 110–114 (2017).
[6] J.-Y. Pan, H.-J. Yang, C. Faloutsos, P. Duygulu, “Automatic multimedia cross-modal correlation discovery” in Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’04 (ACM, New York, 2004), pp 653–658.
[7] L. Cowen, T. Ideker, B. J. Raphael, R. Sharan, Network propagation: A universal am- plifier of genetic associations. Nat. Rev. Genet. 18, 551–562 (2017).