scATAC-seq数据分析之标准分析流程

生信菜鸟团

发布于 2025-05-18 01:11:13

14801

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：1

代码可运行

在上一篇推文【scATAC-seq数据分析之数据读入及质控】中，我们介绍了如何使用 R 包 Signac 进行 scATAC-seq 数据的读取与质控。本节将继续讲解如何基于 Signac 执行标准的分析流程。

首先是加载依赖包：

library(Signac)
library(Seurat)
library(GenomicRanges)
library(ggplot2)
library(patchwork)

Step4. 标准化与线性降维

标准化： Signac 采用 TF-IDF（词频-逆文档频率） 方法进行标准化。这个过程分两步：一是跨细胞标准化，以校正不同细胞测序深度的差异；二是跨峰值标准化，使稀有峰获得更高的权重。
特征选择： 由于 scATAC-seq 数据动态范围较低，因此不像 scRNA-seq 那样进行变异特征选择。可以选择使用前 n% 的特征（峰）用于降维，也可以用 FindTopFeatures() 函数移除出现在少数细胞中的特征。这里我们使用全部特征，但如果只用一部分（比如设置 min.cutoff='q75' 使用前 25% 的峰），运行速度会更快且结果相似。 FindTopFeatures() 会自动设置 Seurat 对象的 VariableFeatures。
降维： 随后我们对 TF-IDF 矩阵进行奇异值分解（SVD），得到降维表示。对于熟悉 scRNA-seq 的用户，可以将其类比为 PCA 降维输出。

TF-IDF 加 SVD 的联合步骤称为 潜在语义索引（LSI），最早由 Cusanovich 等人（2015）用于 scATAC-seq 分析。

pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc)

注意：LSI 第一主成分往往主要反映测序深度（技术变异）而非生物学变异。

可以用 DepthCor() 检查每个 LSI 成分与测序深度的相关性：

DepthCor(pbmc)

image-20250426235012409

在这里，我们观察到第一个LSI主成分与每个细胞的总片段数量之间存在非常强的相关性。因此，在后续分析中，我们将不使用这一成分，因为我们希望根据细胞类型特异性峰的染色质可及性模式来对细胞进行分组，而不是仅仅基于测序深度的差异进行分组。

Step5. 非线性降维与聚类分析

在获得低维表示后，可以像分析 scRNA-seq 一样进行基于图的聚类和非线性降维（如 UMAP）。

pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)

DimPlot(object = pbmc, label = TRUE) + NoLegend() + ggtitle("UMAP of PBMC scATAC-seq data")

Step6. 创建基因活性矩阵

UMAP可视化结果揭示了人类血液中存在多个细胞群体。如果你熟悉外周血单核细胞（PBMC）的scRNA-seq分析，你可能会在scATAC-seq数据中识别出一些髓系细胞群和淋巴系细胞群的存在。

然而，在scATAC-seq数据中对聚类进行注释和解释要更加困难，因为我们对非编码基因组区域的功能角色了解得远远少于对蛋白编码区域（基因）的了解。

为了更好地解析这些细胞群体，我们可以尝试通过评估与基因相关联的染色质开放性，来量化每个基因的活性水平，并基于scATAC-seq数据创建一个新的基因活性矩阵（Gene Activity Assay）。

这里，我们将采用一种简单的方法：统计每个细胞中与基因主体（gene body）及其启动子区域（promoter region）相交的片段数量。（我们也推荐进一步探索 Cicero 工具，它可以实现类似的目标，且在Signac工作流中也提供了运行Cicero的教程。）

为了创建基因活性矩阵，我们需要提取基因的坐标，并向上游延伸2kb区域，因为启动子区域的开放性通常与基因表达水平密切相关。随后，使用 FeatureMatrix() 函数，我们可以统计每个细胞在这些区域内对应的片段数量。

gene.activities <- GeneActivity(pbmc)

pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
  object = pbmc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(pbmc$nCount_RNA)
)

现在，我们可以对一些经典标志基因的活性进行可视化，从而辅助解释ATAC-seq数据中的细胞聚类结果。需要注意的是，相较于scRNA-seq中直接测量的基因表达水平，基于染色质可及性推断的基因活性信号往往更加嘈杂。这是因为基因活性矩阵来源于稀疏的染色质开放性数据，同时也基于一个假设，即基因体及其启动子区域的开放程度能够反映基因的转录活性，但这种对应关系并不总是准确存在。尽管如此，我们依然可以通过基因活性特征，大致辨别出单核细胞、B细胞、T细胞和NK细胞等主要免疫细胞群体。然而，仅依靠这种方式进行细胞亚群的进一步细分仍然具有一定挑战，因此在后续分析中往往需要结合更多信息进行综合判断。

DefaultAssay(pbmc) <- 'RNA'

FeaturePlot(
  object = pbmc,
  features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3
)

image-20250427000134123

Step7. 与 scRNA-seq 数据整合（标签转移）

为了更好地注释 scATAC-seq 数据，可以使用来自相同体系的 scRNA-seq 数据进行标签转移。

加载 10x 提供的预处理 PBMC scRNA-seq 数据：

# Load the pre-processed scRNA-seq data for PBMCs
pbmc_rna <- readRDS("pbmc_10k_v3.rds")
pbmc_rna <- UpdateSeuratObject(pbmc_rna)

通过 跨模态锚点（FindTransferAnchors）寻找对应关系：

transfer.anchors <- FindTransferAnchors(
  reference = pbmc_rna,
  query = pbmc,
  reduction = 'cca'
)

预测 scATAC-seq 细胞的类型：

predicted.labels <- TransferData(
  anchorset = transfer.anchors,
  refdata = pbmc_rna$celltype,
  weight.reduction = pbmc[['lsi']],
  dims = 2:30
)

绘图对比：

pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)
plot1 <- DimPlot(
  object = pbmc_rna,
  group.by = 'celltype',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')

plot2 <- DimPlot(
  object = pbmc,
  group.by = 'predicted.id',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')

plot1 + plot2

image-20250427000222289

处理极少数异常细胞

绘制分配给每个标签的细胞的预测分数揭示了“血小板”细胞得分相对较低（ 0.8），表明对指定细胞身份的信心较低。在大多数情况下，预测这些细胞的下一个最可能的细胞身份是 “CD4 T naive”。如“血小板”这样的极少类群可能是注释误差：

VlnPlot(pbmc, 'prediction.score.max', group.by = 'predicted.id')

image-20250427000243185

# Identify the metadata columns that start with "prediction.score."
metadata_attributes <- colnames(pbmc[[]])
prediction_score_attributes <- grep("^prediction.score.", metadata_attributes, value = TRUE)
prediction_score_attributes <- setdiff(prediction_score_attributes, "prediction.score.max")

# Extract the prediction score attributes for these cells
predicted_platelets <- which(pbmc$predicted.id == "Platelet")
platelet_scores <- pbmc[[]][predicted_platelets, prediction_score_attributes]

# Order the columns by their average values in descending order
ordered_columns <- names(sort(colMeans(platelet_scores, na.rm = TRUE), decreasing = TRUE))
ordered_platelet_scores_df <- platelet_scores[, ordered_columns]

print(ordered_platelet_scores_df)

image-20250427000301379

因此去除总数小于20的预测类型：

predicted_id_counts <- table(pbmc$predicted.id)

# Identify the predicted.id values that have more than 20 cells
major_predicted_ids <- names(predicted_id_counts[predicted_id_counts > 20])
pbmc <- pbmc[, pbmc$predicted.id %in% major_predicted_ids]

# change cell identities to the per-cell predicted labels
Idents(pbmc) <- pbmc$predicted.id

或者，在我们希望结合 scATAC-seq 聚类和标签传输结果的情况下，我们可以将每个 UMAP 聚类内所有细胞的身份重新分配给该聚类最常见的预测标签：

# replace each cluster label with its most likely predicted label
for(i in levels(pbmc)) {
  cells_to_reid <- WhichCells(pbmc, idents = i)
  newid <- names(which.max(table(pbmc$predicted.id[cells_to_reid])))
  Idents(pbmc, cells = cells_to_reid) <- newid
}

Step8. 细胞类型间差异可及性（DA）分析

if (!requireNamespace("remotes", quietly = TRUE))
  install.packages('remotes')
remotes::install_github('immunogenomics/presto')

设置默认 assay 为 peaks，进行差异分析（Wilcoxon秩和检验）：

# change back to working with peaks instead of gene activities
DefaultAssay(pbmc) <- 'peaks'

# wilcox is the default option for test.use
da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  ident.2 = "CD14+ Monocytes",
  test.use = 'wilcox',
  min.pct = 0.1
)

head(da_peaks)
##                           p_val avg_log2FC pct.1 pct.2 p_val_adj
## chr6-13302533-13303459        0  -5.247815 0.026 0.772         0
## chr19-54207815-54208728       0  -4.384378 0.049 0.793         0
## chr17-78198651-78199583       0  -5.505552 0.024 0.760         0
## chr12-119988511-119989430     0   4.153581 0.782 0.089         0
## chr7-142808530-142809435      0   3.663154 0.758 0.087         0
## chr17-82126458-82127377       0   4.997442 0.704 0.043         0

或者，以总碎片数为潜在变量的逻辑回归进行差异表达分析：

另一种鉴别测试的方法是利用逻辑回归，正如 Ntranos 等人所建议的那样 al. 2018 用于 scRNA-seq 数据，并添加片段总数作为潜在变量，以减轻差异测序深度对结果的影响。代码如下：

# change back to working with peaks instead of gene activities
DefaultAssay(pbmc) <- 'peaks'

# Use logistic regression, set total fragment no. as latent var
lr_da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = "CD4 Naive",
  ident.2 = "CD14+ Monocytes",
  test.use = 'LR',
  latent.vars = 'nCount_peaks',
  min.pct = 0.1
)
head(lr_da_peaks)

可视化差异基因：

plot1 <- VlnPlot(
  object = pbmc,
  features = rownames(da_peaks)[1],
  pt.size = 0.1,
  idents = c("CD4 Naive","CD14+ Monocytes")
)
plot2 <- FeaturePlot(
  object = pbmc,
  features = rownames(da_peaks)[1],
  pt.size = 0.1
)

plot1 | plot2

image-20250427000524905

峰值坐标比较难理解。我们可以使用 ClosestFeature() 函数找到与这些峰最接近的基因：

open_cd4naive <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])
open_cd14mono <- rownames(da_peaks[da_peaks$avg_log2FC < -3, ])

closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive)
closest_genes_cd14mono <- ClosestFeature(pbmc, regions = open_cd14mono)

head(closest_genes_cd4naive)
##                           tx_id gene_name         gene_id   gene_biotype type            closest_region              query_region distance
## ENSE00002206071 ENST00000397558    BICDL1 ENSG00000135127 protein_coding exon chr12-119989869-119990297 chr12-119988511-119989430      438
## ENST00000632998 ENST00000632998     PRSS2 ENSG00000275896 protein_coding  utr  chr7-142774509-142774564  chr7-142808530-142809435    33965
## ENST00000665763 ENST00000665763    CCDC57 ENSG00000176155 protein_coding  gap   chr17-82101867-82127691   chr17-82126458-82127377        0
## ENST00000645515 ENST00000645515  ATP6V0A4 ENSG00000105929 protein_coding  cds  chr7-138752625-138752837  chr7-138752283-138753197        0
## ENST00000545320 ENST00000545320       CD6 ENSG00000013725 protein_coding  gap   chr11-60971915-60987907   chr11-60985909-60986801        0
## ENST00000603548 ENST00000603548     FBXW7 ENSG00000109670 protein_coding  utr  chr4-152320544-152322880  chr4-152100248-152101142   219401

head(closest_genes_cd14mono)
##                           tx_id gene_name         gene_id   gene_biotype type           closest_region             query_region distance
## ENST00000606214 ENST00000606214    TBC1D7 ENSG00000145979 protein_coding  gap   chr6-13267836-13305061   chr6-13302533-13303459        0
## ENST00000448962 ENST00000448962      RPS9 ENSG00000170889 protein_coding  gap  chr19-54201610-54231740  chr19-54207815-54208728        0
## ENST00000592988 ENST00000592988     AFMID ENSG00000183077 protein_coding  gap  chr17-78191061-78202498  chr17-78198651-78199583        0
## ENST00000336600 ENST00000336600  C6orf223 ENSG00000181577 protein_coding  utr   chr6-44003127-44007612   chr6-44058439-44059230    50826
## ENSE00002618192 ENST00000569518    VCPKMT ENSG00000100483 protein_coding exon  chr14-50108632-50109713  chr14-50038381-50039286    69345
## ENSE00001389095 ENST00000340607     PTGES ENSG00000148344 protein_coding exon chr9-129752887-129753042 chr9-129776928-129777838    23885

后续我们可以通过对 ClosestFeature() 返回的基因集进行富集分析。

Step9. 差异峰关联基因

我们可以使用 CoveragePlot() 函数，根据细胞的聚类、细胞类型或对象中存储的其他元数据信息，将细胞分组，并绘制Tn5酶切位点在基因组不同区域的分布频率。这种绘图方式生成的是伪bulk染色质可及性轨迹图（pseudo-bulk accessibility tracks），即将同一组中所有细胞的信号进行平均，以可视化特定基因组区域内的DNA开放程度。

除了可及性轨迹之外，我们还可以在同一图中同时展示其他重要信息，例如基因注释、峰的位置信息，以及基因组内的连接关系（如果这些信息已存储在对象中）。更多关于可视化的细节，可以参考官方提供的可视化教程（visualization vignette，https://stuartlab.org/signac/articles/visualization）。

为了使绘图效果更加清晰、易于理解，通常希望将相关细胞类型相邻排列。我们可以通过运行SortIdents() 函数，按照已注释细胞类型之间的相似性，自动调整绘图时的排列顺序。

pbmc <- SortIdents(pbmc)

regions_highlight <- subsetByOverlaps(StringToGRanges(open_cd4naive), LookupGeneCoords(pbmc, "CD4"))
CoveragePlot(pbmc, region = "CD4", region.highlight = regions_highlight, extend.upstream = 1000, extend.downstream = 1000)

image-20250503214715010

regions_highlight <- subsetByOverlaps(StringToGRanges(open_cd14mono), LookupGeneCoords(pbmc, "LYZ"))
CoveragePlot(pbmc, region = "LYZ", region.highlight = regions_highlight, extend.upstream = 1000, extend.downstream = 5000)