scATAC-seq分析ArchR（三）：使用模块打分注释细胞亚群

生信技能树

发布于 2025-05-12 02:53:28

10700

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

scATAC-seq数据分析除了经典的 Signac软件，还有一款也适用差多的软件 ArchR，官网给到你：

https://www.archrproject.com/index.html

这个软件于2021年2月25号发表在 Nature Genetics 杂志上，文献标题为《ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis》。作者认为 ArchR 是目前最全面且用户友好的单细胞ATAC-seq和多组学数据分析工具。开发团队将继续改进和扩展ArchR的功能，以保持其在该领域的领先地位。

大家看文章的通讯作者就知道是 ATAC-seq 这个领域的大佬了：William J. Greenleaf

ArchR分析流程：

续上一篇：

ATAC-seq领域超级大佬William J. Greenleaf团队开发的scATAC-seq分析软件：ArchR

大佬William J. Greenleaf团队开发的scATAC-seq分析软件：ArchR（二）

加载环境

rm(list=ls())
library(ArchR)
set.seed(1)

addArchRThreads(threads = 60) 
## Setting default number of Parallel threads to 1.

addArchRGenome("hg19")
## Setting default genome to Hg19.
# addArchRGenome("hg38")


####### 示例数据下载
library(parallel)
inputFiles <- getTutorialData("Hematopoiesis")
inputFiles
names(inputFiles)

将上一期运行后的数据对象为 projHeme2：经过了降为聚类分群

projHeme2

# class: ArchRProject 
# outputDirectory: /nas2/zhangj/project/12-scATAC-seq/ArchR/Save-ProjHeme1 
# samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
# sampleColData names(1): ArrowFiles
# cellColData names(17): Sample TSSEnrichment ... Clusters ScranClusters
# numberOfCells(1): 10046
# medianTSS(1): 16.8455
# medianFrags(1): 2971

基因打分

ArchR能够稳健地识别聚类，但无法事先确定每个聚类代表的细胞类型，细胞类型的注释通常需要手动完成。现有的自动化注释工具大多为scRNA-seq设计，不完全适用于scATAC-seq数据。

当我们只有 scATAC-seq 数据时，为了进行细胞类型的注释，我们使用对细胞类型特异性标记基因的先验知识，并通过使用基因得分从我们的染色质可及性数据中估计这些基因的表达。基因得分本质上是基于基因周围调控元件的可及性来预测基因的表达水平。

除了基因得分，我们还可以使用模块得分，可以将多个特征（例如基因或峰）组合在一起，以创建更精细的细胞分组定义。

1.计算基因打分

作者测试了50多种不同的基因得分模型，并确定了一类在多种测试条件下始终优于其他模型的模型。这一类模型在ArchR中作为默认模型，它有三个主要组成部分：

基因体内的可及性：对基因得分有贡献；
指数加权函数：以距离依赖的方式考虑潜在远端调控元件的活性；
强制基因边界：最小化不相关调控元件对基因得分的贡献。

在创建Arrow文件时，默认会计算基因得分（addGeneScoreMat参数设置为TRUE），如果在创建时未添加基因得分，可以随时使用addGeneScoreMatrix()函数将基因得分添加到Arrow文件中。基因得分存储在每个Arrow文件中名为“GeneScoreMatrix”的矩阵中。

2.Marker特征鉴定

无偏标记特征识别：ArchR能够无偏地识别任何给定细胞分组的标记特征，这些特征可以是峰、基因（基于基因得分）或转录因子模体（基于chromVAR偏差）。ArchR通过getMarkerFeatures()函数实现这一点，该函数可以通过useMatrix参数接受任何矩阵作为输入，并识别由groupBy参数指示的组所特有的特征。这为无偏地查看每个聚类中预测为活性的基因提供了一种方法，并有助于聚类注释。

标记特征识别的过程如下：

这个函数返回一个SummarizedExperiment对象：

markersGS <- getMarkerFeatures(
  ArchRProj = projHeme2, 
  useMatrix = "GeneScoreMatrix", 
  groupBy = "Clusters",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon"
)

markersGS

# class: SummarizedExperiment 
# dim: 23127 12 
# metadata(2): MatchInfo Params
# assays(7): Log2FC Mean ... AUC MeanBGD
# rownames(23127): 1 2 ... 23126 23127
# rowData names(6): seqnames start ... name idx
# colnames(12): C1 C2 ... C11 C12
# colData names(0):

markerList <- getMarkers(markersGS, cutOff = "FDR <= 0.01 & Log2FC >= 1.25")
head(markerList)
markerList$C6

第6群的特征基因：

可视化：这里挑选了一些已知基因，并将其展示在亚群特征中

markerGenes  <- c(
  "CD34", #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", "EBF1", "MME", #B-Cell Trajectory
"CD14", "CEBPB", "MPO", #Monocytes
"IRF8", 
"CD3D", "CD8A", "TBX21", "IL7R"#TCells
)

heatmapGS <- markerHeatmap(
  seMarker = markersGS, 
  cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
  labelMarkers = markerGenes,
  transpose = TRUE
)
ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")

# 改变cluster顺序：
heatmapGS@row_order <- c(10,11,12,3,4,5,1,2,6,8,7,9)
ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")

# 保存
plotPDF(heatmapGS, name = "GeneScores-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme2, addDOC = FALSE)

3.umap散点图展示基因打分

可以在UMAP嵌入图上叠加每个细胞的基因得分，以便更直观地展示基因表达模式。

markerGenes  <- c(
  "CD34",  #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", "MME", #B-Cell Trajectory
"CD14", "MPO", #Monocytes
"CD3D", "CD8A"#TCells
)

p <- plotEmbedding(
  ArchRProj = projHeme2, 
  colorBy = "GeneScoreMatrix", 
  name = markerGenes, 
  embedding = "UMAP",
  quantCut = c(0.01, 0.95),
  imputeWeights = NULL
)
# 查看单个基因
p$CD14

展示所有基因：

# 展示所有基因
p2 <- lapply(p, function(x){
  x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
      axis.text.x=element_blank(), 
      axis.ticks.x=element_blank(), 
      axis.text.y=element_blank(), 
      axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

# 保存
plotPDF(plotList = p, 
    name = "Plot-UMAP-Marker-Genes-WO-Imputation.pdf", 
    ArchRProj = projHeme2, 
    addDOC = FALSE, width = 5, height = 5)

这种特征基因打分好像并不是很特异啊！

4.MAGIC预测Marker基因打分

单细胞ATAC-seq数据的稀疏性导致基因得分图变化很大。我们可以使用MAGIC通过在邻近细胞之间平滑信号来推算基因得分，这种方法可以改善基因得分的视觉解释。

markerGenes  <- c(
  "CD34",  #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", "MME", #B-Cell Trajectory
"CD14", "MPO", #Monocytes
"CD3D", "CD8A"#TCells
)

p <- plotEmbedding(
  ArchRProj = projHeme2, 
  colorBy = "GeneScoreMatrix", 
  name = markerGenes, 
  embedding = "UMAP",
  imputeWeights = getImputeWeights(projHeme2)
)
p$CD34
p$CD14
p$CD3D

这种方法的基因打分就变得特异了许多：

全部的图并保存：

#Rearrange for grid plotting
p2 <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
do.call(cowplot::plot_grid, c(list(ncol = 3),p2))

plotPDF(plotList = p, 
    name = "Plot-UMAP-Marker-Genes-W-Imputation.pdf", 
    ArchRProj = projHeme2, 
    addDOC = FALSE, width = 5, height = 5)

5.模块打分

从染色质可及性数据中得出的基因得分的可靠性高度依赖于被查询的特定基因，例如，存在于基因组中基因非常密集区域的基因可能更难用基因得分准确表示，因为周围的峰可能并不调节其基因表达。因此，通常最好使用多个不同基因的基因得分来（例如）分配聚类身份。

用addModuleScore()函数添加模块得分，关键参数：

features：一个包含基因名称的命名字符向量列表；
name：将被赋予添加到cellColData的每一列的宽泛标题。

这里，“BScore”模块包含3个B细胞的基因标记（MS4A1、CD79A和CD74），而“TScore”模块包含5个T细胞的基因标记（CD3D、CD8A、GZMB、CCR7和LEF1）。

#############################
features <- list(
  BScore = c("MS4A1", "CD79A", "CD74"),
  TScore = c("CD3D", "CD8A", "GZMB", "CCR7", "LEF1")
)
features

projHeme2 <- addModuleScore(projHeme2,
                            useMatrix = "GeneScoreMatrix",
                            name = "Module",
                            features = features)
                            
# 打分可视化
p1 <- plotEmbedding(projHeme2,
                    embedding = "UMAP",
                    colorBy = "cellColData",
                    name="Module.BScore",
                    imputeWeights = getImputeWeights(projHeme2))

p2 <- plotEmbedding(projHeme2,
                    embedding = "UMAP",
                    colorBy = "cellColData",
                    name="Module.TScore",
                    imputeWeights = getImputeWeights(projHeme2))

ggAlignPlots(p1,p2,draw=T,type="h")

结果如下：

结合umap聚类图，可以很清楚的看到cluster1为B细胞，cluster7,8,9为T细胞：

6.Track Plotting

除了将每个细胞的基因得分作为UMAP叠加图进行绘制外，还可以通过基因组浏览器轨迹在每个聚类的基础上浏览标记基因的局部染色质可及性。使用plotBrowserTrack()函数，该函数将为每个指定的基因创建一个绘图列表：

########################################
set.seed(1)

markerGenes  <- c(
"CD34", #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", #B-Cell Trajectory
"CD14", #Monocytes
"CD3D", "CD8A", "TBX21", "IL7R"#TCells
)

p <- plotBrowserTrack(
  ArchRProj = projHeme2, 
  groupBy = "Clusters", 
  geneSymbol = markerGenes, 
  upstream = 50000,
  downstream = 50000
)
p
grid::grid.newpage()
grid::grid.draw(p$CD14)