scATAC揭示了单细胞分辨率的表观遗传景观。 该测定法已成功应用于鉴定细胞类型及其特定的调控件揭示细胞异质性，映射疾病相关的远端元素并重建分化轨迹 。然而，由于accessible peaks的固有高维度，以及单细胞sequencing reads的稀疏性，分析scATAC 数据仍然面临着巨大的挑战。

概述

scATAC数据具有以下几个特点：

  • 高维度

    检测到的可能开放区域高达几十万

  • Binary

    每个开放区域在二倍体基因组上通常只有两个拷贝，导致数据接近二值化

  • 极端稀疏性

    由于细胞异质性以及技术原因，导致每个单细胞都有大量的开放区域没有信号

传统的降维和可视化方法以及基于样本距离聚类的方法无法在单细胞ATAC-seq数据上取得很好的效果。且由于单细胞ATAC-seq相比于单细胞RNA-seq数据，具有二值化和更加稀疏的特性，直接应用单细胞RNA-seq的方法也不太理想。

分析流程：

  • 质控 -- 参考ENCODE
    • TF可及性分析
    • 共可及性分析
    • GeneScore Matrix
  • 聚类
    • 策略1，Tile Matrix
    • 策略2，PeakMatrix
    • 策略3，TF-motif Matrix
  • 细胞类型注释
    • 策略1，Marker Gene可及性人工注释
    • 策略2，GeneScore人工/软件注释
    • 策略3，类间特异的TF-motif
    • 策略4，整合scRNA
    • 筛选细胞类型
      • TF/Peak可及性差异分析
      • 共可及性分析
      • 拟时间分析

质控

fragment

如cellranger-atac算法中，认为fragments数量超过阈值，且保守可信的peaks区域中的fragments的比例大于阈值时，barcodes才是有效细胞。

TSS enrichment score

基因转录时开放TSS区域，结合转录因子等转录起始复合物，一般情况下TSS区域的fragments占细胞全部fragments的25%-35%。一般要求有效细胞barcode的TSS富集得分>4，且barcode包含fragments>1000。

相当于计算每个细胞的信噪比：基于TSS中心的fragments和TSS侧翼区域fragments的比例。

nucleosome_signal

DNA开放区域长度，一般在1-2个核小体范围内，一般过滤掉nuclesome-signal>4的barcodes

blacklist regions ratio

blacklist regions 基因组中总反常或二代测序中的高信号区域。过滤掉

双细胞

一个barcode中包含的fragments太多，说明可能包含两个或多个细胞。fragments数目阈值或demuxlet双细胞过滤，对Doublet Enrichment得分从大到小排序，按比例过滤。

聚类

策略1，Tile Matrix

设置滑动窗口bin（如5kb），依次统计单个细胞在每个bin中fragments的数量。

fragments指的是atac-seq双端测序后得到的序列，分析过程中的fragment指的是记录在染色体上，经过偏移校正后，单碱基的起始位置和结束位置。同样insertion指的是偏移校正后开放位置中心的单碱基位置。

策略2，PeakMatrix

对fragment进行peak calling，过滤可信度低的数据（在少数细胞中出现的peak），最终生成Peak-Barcode可及性矩阵。在marker gene上下游特别是TSS区域Peaks分布差异，知道细胞分类。scRNA数据feature数目大约2-3万，scATAC peaks大约10万左右。

基于上述矩阵进行降维

策略3，TF-motif Matrix

一个TF通常包含一个可结合的DNA motif，少数情况下有两个或多个。对peaks反向搜索，统计每个TF-motif在每个细胞中可以结合的fragments数目，得到单细胞层面、转录因子水平的可及性矩阵TF-motif Matrix。chromeVAR。 TF-motif库，中数据在600左右，所以从Peak-motif可以看作数据降维的过程。

还有基于Gene Score Matrix进行聚类

初始聚类->peak calling->再聚类

  • 有效细胞筛选

  • TileMatrix聚类

  • 对cluster中，每200个细胞合并为一个pseudo bulk ATAC数据

  • 对每个pseudo bulk ATAC数据做call peaks，再合并

  • 得到Peak Matrix

合并为pseudo bulk ATAC数据，有效减少scATAC数据捕获不足导致的数据缺失；同时callpeak得到细胞群间特异的peaks。

peaks标准化

log（TF-IDF）文档频率法：TF是词频(Term Frequency)，IDF是逆文本频率指数(Inverse Document Frequency)。单个词汇在一篇文章中出现的次数越多，越重要。但是在语料库多次出现，重要性越来越低。IDF : 计算A term 出现稀少度。越稀少，越重要。

peaks降维

LSI：潜在语义索引

细胞类型注释

Maker基因

在细胞分群结果基础上，进行类间peak可及性分析，根据marker基因TSS区域或整个基因区域的可及性类间差异，对细胞群进行人工注释。

基因活性

Gene Score Matrix基因活性得分矩阵，考虑peak与TSS/gene body的距离、测序深度等因素采用广义线性模型将Peak Matrix转换为Gene Score Matrix，常用软件Cicero、Signac、ArchR。

基因活性得分可以理解为基于表观数据预测出来的基因表达值。可以通过观察marker gene在细胞群中活性得分值的特异分布，进行细胞类型人工注释、或SingleR、Cellassign

联合scRNA转移类标签

寻找最接近且已经注释好细胞类型的scRNA数据，基于Gene Score，采用Seurat TransferData方法，获得scATAC的细胞类型注释

共可及性分析

peak-peak共可及性分析，常用来研究基因组上的顺式作用元件间的相互作用，Cicero

横坐标预测到的共可及性网络区间内的peaks和注释到的基因，纵坐标共可及性值，弧线为peak之间的相互作用连接线

拟时间分析（趋势分析）

基于基因活性得分矩阵，Monocle软件实现scATAC数据的拟时间分析

一篇单细胞ATAC-seq分析的综述：
Shi, P., Nie, Y., Yang, J. et al. Fundamental and practical approaches for single-cell ATAC-seq analysis. aBIOTECH 3, 212–223 (2022). https://doi.org/10.1007/s42994-022-00082-5