单细胞分析之细胞交互-3：CellChat

Posted 2023-05-04

参考技术A

CellChat 通过综合信号配体、受体及其辅因子基因的表达只与它们之间互作的先验知识对细胞通讯概率建模。在推断出细胞间通讯网络后，CellChat提供了进一步数据探索、分析和可视化的功能。

文章链接： https://www.nature.com/articles/s41467-021-21246-9
视频链接： https://www.youtube.com/watch?v=kc45au1RhNs

CellChat工作流程图：

从Seurat对象直接创建：
⚠️：构建Cell Chat对象时，输入的是log后的数据。

在CellChat中，我们还可以先择特定的信息描述细胞间的相互作用，可以理解为从特定的侧面来刻画细胞间相互作用，比用一个大的配体库又精细了许多。

对表达数据进行预处理，用于细胞间的通信分析。首先在一个细胞组中识别过表达的配体或受体，然后将基因表达数据投射到蛋白-蛋白相互作用(PPI)网络上。如果配体或受体过表达，则识别过表达配体和受体之间的相互作用。

通过为每个相互作用分配一个概率值并进行置换检验来推断生物意义上的细胞-细胞通信。

通过计算与每个信号通路相关的所有配体-受体相互作用的通信概率来推断信号通路水平上的通信概率。

我们可以通过计算链路的数量或汇总通信概率来计算细胞间的聚合通信网络。

至此，统计分析部分已经完成。

检查每种细胞发出的信号

⚠️注：层次图、网络图、和弦图、热图只是不同的展示方法，展示的内容和代表的意思一模一样
比如在前面的功能富集分析或case control的比较中找到了一些信号通路差异，就可以进一步在细胞水平上验证。

层次图（Hierarchy plot）

网络图（Circle plot）

和弦图（Chord diagram）

热图（Heatmap）

层次图（ Hierarchy plot）

网络图（Circle plot）

和弦图（Chord diagram）

Error: not enough space for cells at track index \'1\'.

气泡图（全部配体受体）

气泡图（指定信号通路或配体-受体）

气泡图（指定信号通路或配体-受体并指定细胞）

参与某条信号通路（如TGFb）的所有基因在细胞群中的表达情况展示（小提琴图和气泡图）

通讯网络系统分析使用了三种算法： 社会网络分析 、 NMF分析 和 流行学习与分类
⚠️：不同的算法算出来的结果可能会相互矛盾，需要结合生物学知识加以判断

计算网络中心性权重

通过计算每个细胞群的网络中心性指标，识别每类细胞在信号通路中的角色/作用C（发送者、接收者、调解者和影响者）

识别细胞的信号流模式

热图

冲击图/河流图（river plot）

气泡图

热图

冲击图

气泡图

把共同起作用的信号通路归纳在一起，分为基于功能的归纳和基于拓扑结构的归纳

⚠️：配对分析必须保证细胞类型是一样的，才可以进行配对。如果 两个样本的细胞类型不一样又想进行配对分析时，可以用subset把两个样本的细胞类型取成一致的。

数量与强度差异网络图

数量与强度差异热图

细胞互作数量对比网络图

这里function的图出不来，只有structural的图可以出来

总体信号模式对比

单细胞ATAC亚群分析

参考技术A

上一期，跟大家简单介绍了 下单细胞ATAC的背景知识点及其10x ATAC基础数据的获取方式 。接下来就带大家从fragment.csv、singlecell.csv、peaks matrix等数据出发，做单细胞ATAC的亚群分析。

与单细胞转录组类似，单细胞ATAC的分析流程也主要包括细胞质控、peaks标准化及其降维分群、marker基因的鉴定等几个步骤。常用的单细胞ATAC分析流程软件包含 cell-ranger-atac、Signac和ArchR等。

单细胞ATAC的质控点一般包含以下几个方面：样本重复（biological replicates），bulkATAC vs scATAC的相关性、fragment length distribution、per nucleus read-depth、transcription start site (TSS) enrichment、双细胞比例等。

前面提到的样本相关性和fragments的长度分布主要是从整体水平上检查我们的单个样本数据的可靠性。
而要去掉不符合质控的细胞，我们主要从 fragments 数目 和 TSS enrichment score 这两点出发。

👉 fragments 数目：一般指单个细胞（barcode）所属的total fragments数目。这个不同的软件具体的定义不同，比如cell-ranger-atac和Signac指peaks所属区域的fragments 数目，其中singlecell.csv文件中peak_region_fragments列便是指fragments 数目，而ArchR是指全基因组所有的fragments 数（这个跟该软件的分析策略有关，后面会提到）。

备注：fragments 数目&TSS enrichment score的阈值不仅与所用软件具体的计算公式有关（不同的软件具体的参数可能不同），也与自己数据的实际情况有关。比如哺乳动物和植物的单细胞ATAC数据TSS enrichment score就不能用相同的指标cutoff来衡量，一般来说哺乳动物的TSS enrichment score值要整体偏高些。

双细胞预测几乎是所有单细胞测序技术都得考虑的一个问题，从原理来说，我们每个barcode就是一个细胞，但是因为所有的实验技术都不是100%完美的，因此往往会有一个barcode所包裹的油滴进来2个细胞。
对于10x数据来说，即使在使用标准试剂盒时，也可能有超过5%的细胞属于双细胞，这对聚类产生了重大影响。特别是在发育/轨迹分析中十分受影响，因为doublets看起来像是两种细胞类型的混合物，这可能与中间细胞类型或细胞状态混淆。

为了预测哪些“细胞”实际上是双细胞的，ArchR会从我们真实的数据中随机模拟产生混合的“双细胞”数据，这些“双细胞”数据与我们所有细胞一起做降维并UMAP可视化（"双细胞"会投影到UMAP中，并识别它们邻近的细胞），在这个过程中，ArchR会计算每个细胞的Doublet Enrichment，值越大，表示该细胞是双细胞的可能性越大。

与单细胞RNA(scRNA-seq)相比，scATAC-seq数据由于其高维度和稀疏性而更具计算分析挑战性。主要体现在标准化和降维，这两大步骤跟单细胞转录组分析所用的统计学原理完全不同，以下为归纳总结的具体内容，如下表所示：

备注：TF-IDF & LSI都是自然语言常用的统计学方法。

获得peak matrix后，跟基因类似，我们必须对其标准化。因为单细胞ATAC测的是DNA序列，对于二倍体物种来说，同一个位置最多有2套DNA序列，这便是单细胞ATAC peak matrix稀疏性的最大根源（单细胞转录组因测的是RNA，高表达的基因往往有多个转录分子）。因此，从数据实际情况出发，单细胞ATAC采取的是log(TF-IDF)( Term frequency-inverse document frequency) 标准化，简称文档频率法。

所有高维数据的分析都是采取降维的方式从多维到低纬的策略，之后还可以再次降维成2个维度并可视化（比如TSNE和UMAP）。我们对peaks是采取LSI降维的方式。

与单细胞转录组类似，降维后的单细胞ATAC数据也同样可以采取graph-based clustering的分群方法。Graph-based图聚类算法包括两步：首先用降维（PCA或者LSI）的数据构建一个细胞间的k近邻稀疏矩阵，即将一个细胞与其欧式距离上最近的k个细胞聚为一类，然后在此基础上用Louvain算法进行模块优化(Blondel, Guillaume, Lambiotte, & Lefebvre, 2008)，旨在找到图中高度连接的模块。最后通过层次聚类将位于同一区域内没有差异表达基因(B-H adjusted p-value 低于0.05)的cluster进一步融合，重复该过程直到没有clusters可以合并。
备注：Signac和ArchR都是直接调用Seurat包的FindClusters()函数用不同分辨率来分群的。

细胞分群后，我们需要知道每个cluster属于什么细胞类型，也就是细胞命名。我们知道，单细胞转录组主要是依据每个cluster的marker基因来判断细胞类型的。那么 对于单细胞ATAC，是不是也可以定义出每个cluster的特异高表达的基因集呢？
答案是肯定的，一般来说，我们是通过基因body区域加上一定范围内的上下游区域的整体ATAC信号来计算每个细胞每个基因的genescore值。

1）Signac是通过 GeneActivity() 函数 https://satijalab.org/signac/reference/geneactivit 来实现的，默认参数是基因上游2kb到TES区域。
2）而ArchR是通过 addGeneScoreMatrix() 函数 https://www.archrproject.com/reference/addGeneScoreMatrix.html 来实现的（createArrowFiles函数也会用默认参数得到genescore matrix矩阵），注意其计算原理稍微复杂，ArchR考虑到远端调控元件对基因活性的影响，因此默认的upstream和downstream范围更广。

在ArchR作者的发表文章中，他们测试了50多个不同的基因评分模型，并确定了一类在各种测试条件下表现始终优于其他模型的模型。这个模型类，在ArchR中作为默认实现，有三个主要组件:

marker 基因的ATAC信号（genescore值）同样可以在umap展示，也可以用小提琴图（VlnPlot），点状图（DotPlot）展示。与单细胞转录组相比，单细胞ATAC还多了基因区域的track的可视化展示。

单细胞ATAC的亚群分析介绍就到这里，下一篇会给大家介绍单细胞ATAC的高级分析内容，比如motifdeviation、 拟时间分析、 单细胞RNA与单细胞ATAC的整合分析等。
本分享更多是从知识点和分析原理来讲解和归纳总结，具体实现方法和流程脚本可以查看下面参考资料软件的官方文档，里面都写得都很详细清楚。