scATAC-seq 基因活力分析：染色质可及性与转录潜力评估
时长: 30 分钟
字数: 7.8k 字
更新: 2026-03-02
阅读: 0 次
背景介绍
基因活力分析（Gene Activity Analysis）是单细胞 ATAC-seq 数据分析中的核心技术，它通过量化染色质可及性来间接评估基因的转录潜力。这种方法基于一个重要的生物学假设：染色质的开放状态与基因转录活性之间存在密切关联。
核心原理
该方法通过统计每个基因区域（包括基因体及其上游 2kb 调控区域）内的 ATAC-seq 信号强度，这些信号直接反映了 DNA 的开放程度，从而预示着基因的潜在转录活性。
目的与意义
- 数据整合桥梁：为 scATAC-seq 构建“转录近似”表型，用作与 scRNA-seq 的锚点，支持 CCA/Anchors 等跨组学整合与标签转移。
- 注释与解释：在 ATAC 聚类上可视化经典标记基因活力，辅助识别细胞类型与状态。
- 机制与发现：结合表达与活力做相关/差异分析，挖掘特异开放基因与潜在调控逻辑。
本教程涵盖内容
- 基因活力分析：评估染色质可及性与基因转录潜力的关系
- 基因活力差异分析：识别不同细胞类型间的基因活力差异
- 差异分析可视化：结合基因表达数据，对比差异基因活力与对应基因的表达模式
- 相关性分析：
- 不同 cluster 间的基因活力相关性分析
- 不同 cluster 间的基因表达相关性分析
预期分析结果
通过本部分的分析，您将获得：
- cluster 特异性基因集：识别在特定细胞类型中特异性开放的基因
- 表达-活力关联性：了解 cluster 特异性基因在基因活力和基因表达两个层面的表现
- 异质性评估：评估不同 cluster 在基因表达层面和基因活力层面的异质性特征
参数设置
请先挂载将要分析的云平台相关流程数据，挂载方式请参照 jupyter 使用教程
具体参数填写：
- rds：挂载的流程相关的 rds 数据，在/home/{UserName}/workspace/project/{UserName}/目录下，如下列项目数据/home/shumeng/workspace/data/AY1752167131383/
- meta：与 rds 同目录下的 metadata 文件
- outdir: 分析结果报错路径
- celltype_colname：细胞类型列名，选择要分析的细胞类型或者聚类对应的标签。假如想对注释好的细胞类型进行分析，则选择其对应的标签，例如 CellAnnotation，与<细胞类型>配合使用
- celltypes：细胞类型，多选，选择要分析的细胞类型或者聚类结果，如 Monocyte 和 Macrophage 等
- logFC_cut：差异倍数的阈值，默认 0.25
- padj_cut：pvadj 阈值，默认 < 0.05
- pval_cut：pval 阈值，默认 < 0.05
- reduction：FeaturePlot 可视化基因活力时，用于降维的方式，可选 umap/tsne/lsi/atacumap/atactsne/wnnumap/wnntsne。默认 wnnumap
R
# Parameters按照需要填写自己路径下的数据路径
rds = "/home/shumeng/workspace/data/AY1752167131383/input.rds"
meta = "/home/shumeng/workspace/data/AY1752167131383/meta.tsv"
outdir = "/home/shumeng/workspace/project/shumeng/"
celltype_colname = "wknn_res.0.5_d30_l2_50"
celltypes = "0,1,2,3,4,5,6,7,8,9"
logFC_cut = 0.25
padj_cut = 0.05
pval_cut = 0.05
reduction = "wnnumap"环境准备
R 包加载
请选择 common_r 这个环境进行该教程的学习
R
# 保存当前库路径
original_lib_paths <- .libPaths()
# 1. 首先加载基础依赖包
suppressMessages({
library(Rcpp, quietly = TRUE)
library(Matrix, quietly = TRUE)
})
# 2. 设置临时库路径并加载 presto
temp_lib_path <- "/PROJ/development/liuxin/Software/miniconda3/envs/Celltype/lib/R/library"
if(dir.exists(temp_lib_path)) {
.libPaths(c(temp_lib_path, original_lib_paths))
}
suppressMessages({
library(presto, quietly = TRUE)
source("/PROJ2/FLOAT/shumeng/script/wilcoxauc.R")
})
# 3. 恢复原始库路径
.libPaths(original_lib_paths)
# 4. 加载其他分析包，调整加载顺序
suppressMessages({
# 先加载基础依赖
#library(parallelly, quietly = TRUE)
#library(future, quietly = TRUE)
# 然后加载主要分析包
library(ggrepel,quietly = TRUE)
library(Signac,quietly = TRUE)
library(Seurat,quietly = TRUE)
# 最后加载可视化和工具包
library(ggplot2, quietly = TRUE)
library(dplyr, quietly = TRUE)
#library(ggrepel, quietly = TRUE)
library(DT, quietly = TRUE)
library(hexbin, quietly = TRUE)
library(viridis, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(pheatmap, quietly = TRUE)
library(base64enc, quietly = TRUE)
library(viridis, quietly = TRUE)
library(RColorBrewer,quietly = TRUE)
library(scales,quietly = TRUE)
})基因活力（Gene Activity）计算
计算流程
- 目标区域确定：为每个基因定义基因体（gene body）与上游约 2 kb 启动子区域（可按需要调整）。
- 片段计数：对每个细胞统计 ATAC 片段落在上述区域内的次数（可按片段或插入位点计数）。
- 构建矩阵：将基因×细胞的计数整合为基因活力矩阵（
FeatureMatrix；在实际流程中由GeneActivity()封装完成）。 - 归一化：对活力矩阵进行对数标准化，使其量级可与 scRNA-seq 表达相比较，便于联合分析与可视化。
在本项目中的实现（Signac）
- 读取单样本或整合后的
Seurat对象（可选合并meta并AddMetaData）。 - 设置默认
assay为ATAC，运行GeneActivity(sc)计算基因活力。 - 将结果写入新的
assay：ATAC_activity（CreateAssayObject(counts = gene.activity)）。 - 使用
NormalizeData对ATAC_activity进行LogNormalize，scale.factor取median(sc$nCount_RNA)，以便与 RNA 表达量处于相近量级，支持后续整合（如 CCA/Anchors）与对比展示。
实践要点与注意事项
- 间接性与噪声：活力分数来自稀疏的可及性数据，较表达更噪；与真实表达并非一一对应。
- 区域定义敏感：启动子窗口（如 2 kb）及是否包含基因体会影响结果，可按物种与注释版本微调。
- 远端调控影响：对由远端增强子主导的基因，活力与表达的一致性可能较弱；可结合
Cicero/ArchR等工具建联系并验证。
R
suppressWarnings({
suppressMessages({
sc=readRDS(rds)
if(!is.null(meta)){
meta <- read.table(meta, header = T, sep = "\t")
rownames(meta) <- meta$barcodes
sc <- AddMetaData(sc, meta)
}
DefaultAssay(sc)="ATAC"
gene.activity <- GeneActivity(sc)
sc[['ATAC_activity']] <- CreateAssayObject(counts = gene.activity)
sc <- NormalizeData(
object = sc,
assay = 'ATAC_activity',
normalization.method = 'LogNormalize',
scale.factor=median(as.numeric(sc$nCount_RNA))
)
})
})R
sc@assays$ATAC_activity$counts[1:6,1:6]基因活力矩阵与基因表达矩阵类似，每一行代表一个基因，每一列代表一个细胞。值则表示每个细胞中对应基因的基因活力值。值越大表示该基因在对应细胞中的调控区域越开放，对应基因的潜在转录活性越高。此结论建立在基因活力与基因表达存在正相关的假设之上，实际由于基因表达调控的复杂性，可能存在基因活力高但基因表达低的情况，具体调控关系可根据实际需求进一步探究。基因活力差异分析
细胞分群间基因活力差异分析
分析目的
识别不同细胞分群（cluster）中基因活力显著升高或降低的基因，识别在特定细胞类型中特异性开放的基因，刻画各分群的表观转录特征。
对比策略与方法
- 对比方式：采用“一对其余”（one-vs-rest）策略，对每个分群与其他所有分群进行独立比较。
- 统计检验：使用 Wilcoxon 秩和检验评估基因活力在目标分群与其它分群之间的差异。
输出指标说明
- feature：基因名称（基因层面的活力特征）
- group：目标细胞分群
- avgExpr：目标分群内该基因的平均活力水平
- logFC：log2 倍数变化（目标分群相对其余分群的差异；正值表示目标分群更高）
- statistic：Wilcoxon 检验统计量
- auc：二分类 AUC，表示该基因活力对目标分群的区分能力
- pval：原始 P 值
- padj：多重检验校正后的 P 值
- pct_in / pct_out：该基因在目标分群/其它分群中被检测到（活力>0）的细胞比例
筛选与阈值
可结合效应大小与显著性设定阈值（如对 logFC 与 padj 同时约束），并参考检测率（pct_in/pct_out）以提升结果的稳健性。筛选后可获取各分群的候选高活力（或低活力）基因集。
R
dir.create("../result", recursive = TRUE, showWarnings = FALSE)
dir.create("../result/cluster_Diff", recursive = TRUE, showWarnings = FALSE)
dir.create("../result/cluster_Diff/FeaturePlot", recursive = TRUE, showWarnings = FALSE)
#dir.create("../result/group_Diff", recursive = TRUE, showWarnings = FALSE)R
#celltype或者和cluster之间进行基因活力差异分析
DefaultAssay(sc) <- "ATAC_activity"
suppressWarnings({
suppressMessages({
Idents(sc)=sc@meta.data[[celltype_colname]]
Diff_genes_activity=wilcoxauc(sc,celltype_colname,assay = 'data',seurat_assay = 'ATAC_activity')
if(exists("padj_cut")){
Diff_genes_activity=Diff_genes_activity[which(abs(Diff_genes_activity$logFC) > logFC_cut & Diff_genes_activity$padj < padj_cut),]
}else if(exists("pval_cut")){
Diff_genes_activity=Diff_genes_activity[which(abs(Diff_genes_activity$logFC) > logFC_cut & Diff_genes_activity$padj < pval_cut),]
}
write.table(Diff_genes_activity,file = "../result/cluster_Diff/cluster_Diff_genes_activity.txt",quote = F,sep = "\t",row.names = F)
})
})结果呈现与使用
结果以交互式表格展示，支持过滤、排序与导出，便于提取各分群特异的高活力基因集
R
# 移除外层抑制函数，保证输出可见性
if(nrow(Diff_genes_activity) < 1) {
print("没有基因活力显著差异的基因，请降低显著性阈值")
} else {
# 直接返回DT对象（无需print）
datatable(
Diff_genes_activity,
rownames = FALSE,
filter = "top",
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
buttons = c('csv', 'excel'),
columnDefs = list(
list(
targets = 2:9, # 更简洁的索引写法
render = DT::JS("function(data) {
return typeof data === 'number'
? data.toExponential(2)
: data;
}")
)
)
)
)
}通常当同时满足 logFC > 0.25 和 padj < 0.05 这两个条件时，可以认为该基因在这个特定细胞群中相比其他细胞群具有显著更高的调控活力。降维图可视化差异结果
特征基因可视化（FeaturePlot）
为直观展示不同细胞类型的标志基因分布，我们为每个细胞类型选取最具代表性的前五个标志基因，并在降维空间中进行并行可视化。
- 展示维度
- ATAC_activity：反映基因启动子/基因体区域的染色质开放程度（调控潜力）
- RNA：反映对应基因的转录表达水平（实际转录产物）
- 图形说明
- 每个点代表一个单细胞；颜色深浅表示该基因在该细胞的活力或表达强度
- 将同一基因在两种模态下并列/分面展示，便于对比不同细胞类型中的特异性模式
- 解读要点
- 细胞类型特异性：标志基因在预期细胞类型中应呈现高活力/高表达的空间富集
- 跨模态一致性：ATAC 活力与 RNA 表达的一致性支持该基因在该类型的活跃转录状态
- 不一致线索：若活力高但表达低，可能提示调控预激活、时序差异或远端增强子影响
- 注意事项
- 基因活力为间接读数，噪声较大；与 RNA 表达不一定一一对应
- 活力的窗口定义（如上游 2 kb 是否包含、是否计入基因体）会影响可视化结果
R
#获取每个cluster中top10差异基因，和热图基因，为后续可视化做准备
if(nrow(Diff_genes_activity) < 1){
print("没有基因活力显著差异的基因，请降低显著性阈值")
}else{
cluster_num <- length(table(sc@meta.data[[celltype_colname]]))
if (cluster_num <= 8) {
n <- 10
height <- 3
} else if (cluster_num > 8 & cluster_num <= 10) {
n <- 8
height <- 4
} else if (cluster_num > 10 & cluster_num <= 15) {
n <- 5
height <- 5
} else {
n <- 4
height <- 7
}
Diff_genes_activity %>%
group_by(group) %>%
arrange(desc(logFC)) %>%
slice_head(n = n) %>% # 改为10个基因
ungroup() -> top10
if(nrow(top10) <= 10){
width=5
}else if(nrow(top10) > 10 & nrow(top10) <= 15){
width=7
}else if(nrow(top10) > 15 & nrow(top10) <= 20){
width=9
}else if(nrow(top10) > 20 & nrow(top10) <= 25){
width=11
}else if(nrow(top10) > 25 & nrow(top10) <= 30){
width=13
}else if(nrow(top10) > 30 & nrow(top10) <= 35){
width=15
}else if(nrow(top10) > 35 & nrow(top10) <= 40){
width=17
}else if(nrow(top10) > 40 & nrow(top10) <= 45){
width=19
}else if(nrow(top10) > 45 & nrow(top10) <= 55){
width=20
}else if(nrow(top10) > 55 & nrow(top10) <= 65){
width=21
}else {
width=22
}
#获取每个cluster中top1差异基因
Diff_genes_activity %>%
group_by(group) %>%
arrange(desc(logFC)) %>% # 添加这行，按照logFC从大到小排序
slice_head(n = 5) %>% # 改为1个基因
ungroup() -> top5
#获取画热图的基因集
Diff_genes_activity %>%
group_by(group) %>%
arrange(desc(logFC)) -> heatmap_genes
}R
# 添加必要的包
# 检查并创建保存目录
# dir.create("_build/html/result/cluster_Diff/FeaturePlot/", recursive = TRUE, showWarnings = FALSE)
options(warn = -1)
save_and_encode_plot <- function(plot, filename, width, height) {
# 保存图片文件
full_path_pdf <- paste0(filename, ".pdf")
temp_png <- tempfile(fileext = ".png")
# 保存PDF
ggsave(full_path_pdf, plot, width = width, height = height)
# 临时保存PNG用于base64编码
ggsave(temp_png, plot, width = width, height = height, dpi = 300)
# 将图片转换为base64
png_data <- readBin(temp_png, "raw", file.info(temp_png)$size)
base64_img <- base64encode(png_data)
# 删除临时PNG文件
#unlink(temp_png)
return(list(
base64 = base64_img,
pdf_path = full_path_pdf
))
}
suppressWarnings({
suppressMessages({
if(exists("top5")){
cell_types <- unique(top5$group)
# 创建数据框用于展示
image_df <- data.frame(
CellType = cell_types,
stringsAsFactors = FALSE
)
# 为每个assay生成一列图片
for(assay in c("ATAC_activity", "RNA")) {
DefaultAssay(sc) <- assay
col_name <- paste0(assay, "_Plot")
image_df[[col_name]] <- sapply(cell_types, function(i) {
# 1. 生成图片
p1 <- FeaturePlot(
object = sc,
features = top5$feature[top5$group == i],
reduction = reduction,
pt.size = 0.1,
max.cutoff = 'q95',
ncol = 5,
cols=c("#E5E3E3","#EF6067"),
order = TRUE
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.text = element_text(size = 8),
axis.title = element_text(size = 10)
)
# 2. 保存图片并获取base64编码
result <- save_and_encode_plot(
plot = p1,
filename = paste0("../result/cluster_Diff/FeaturePlot/feature_plot_cluster_",
assay, "_", i),
width = 20,
height = 3
)
# 3. 创建HTML，包含base64图片
sprintf(
'<div class="plot-container">
<img src="/placeholder.svg" height="400px" onclick="window.open(\'%s\')"/>
</div>',
result$base64,
result$pdf_path
)
})
}
# 4. 创建交互式表格展示结果
output_table <- suppressWarnings(DT::datatable(
image_df,
escape = FALSE,
options = list(
pageLength = 1,
scrollX = TRUE,
dom = 'Bfrtip',
server = TRUE, # 强制服务器端模式
pageLength = 1,
deferRender = TRUE, # 延迟渲染
scroller = TRUE, # 滚动加载
columnDefs = list(
list(
targets = 1:2,
width = "45%"
)
)
),
rownames = FALSE,
filter = 'top',
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: center; font-size: 1.2em; font-weight: bold;',
'Feature plots for top 5 genes in each celltype'
)
)
)
output_table
}else{
print("没有基因活力显著差异的基因，请降低显著性阈值")
}
})
})
options(warn = 0)表格中第一列为细胞分群信息，第二列为该细胞群中最具代表性的前五个标志基因染色质活力分布图，第三列为前五个标志基因基因表达分布图热图可视化差异结果
特征基因热图对比（ATAC_activity vs RNA）
可视化目的
- 展示各细胞类型标志基因在两种模态（染色质活性与转录表达）下的分布与一致性，辅助细胞类型注释与机制解读。
图像构成
- 左侧：ATAC 活性热图（ATAC_activity，反映启动子/基因体的可及性）
- 右侧：RNA 表达热图（反映实际转录水平）
- 行：特征基因（每个细胞类型选取前五个或自定义列表）
- 列：单细胞；按细胞类型/cluster 分组，并以颜色条标注
- 颜色：由浅到深表示活性/表达由低到高
- 标准化建议：为突出相对模式，建议对每个基因按细胞方向做行内标准化（如 z-score），并在两张热图中保持相同的基因顺序与列顺序
解读要点
- 细胞类型特异性：标志基因在对应细胞类型中应显示高活性/高表达的聚集带
- 跨模态一致性：ATAC 活性与 RNA 表达的协同上调，提示活跃的转录状态与开放染色质的一致
- 不一致线索：活性高但表达低，可能反映转录预激活、时间滞后或远端增强子驱动；表达高但活性低，可能受技术噪声或区域定义影响
- 异质性：同一分群内的梯度/多峰提示亚群或连续谱状态
R
# 检查并处理 color_gradient
# 将逗号分隔的字符串拆分成颜色值向量
if(!exists("color_gradient")){
color_gradient <- c("#F9F9F9","#F0484A")
} else {
# 如果 color_gradient 是字符串，将其分割成向量
if(length(color_gradient) == 1) {
color_gradient. <- unlist(strsplit(color_gradient, ","))
}
}
colors_cancer <- colorRampPalette(color_gradient)(100)
base_colors=c("#EDD496","#65D1BF","#7ADAEA","#F2AFAF","#EA8F8F",
"#C1A9D3","#F2DA9E","#D9EF78","#EF9ED9","#C9BCC6","#8BD198")
library(RColorBrewer)
library(scales)
# 获取实际的cluster数量
cell_types <- if(exists("celltypes")) unlist(strsplit(celltypes, ",")) else unique(sc@meta.data[[celltype_colname]])
unique_clusters <- length(cell_types)
cell_types_all=cell_types
generate_cols <- function(n) {
if(n <= 0) return(character(0)) # 处理0长度情况
# 循环扩展基础颜色
base_colors[((1:n)) %% length(base_colors) + 1]
}
my_cols <- generate_cols(unique_clusters)
identity_colors <- setNames(my_cols, cell_types)R
suppressWarnings({
suppressMessages({
# 加载必要的包
# 设置图形参数
options(repr.plot.width = 16, repr.plot.height = 10)
if(exists("top10")){
# ATAC活性热图
DefaultAssay(sc) <- 'ATAC_activity'
p1 = DoHeatmap(sc,
features = top10$feature,
slot = "data",
draw.lines = TRUE,
raster = TRUE,
group.colors = identity_colors,
label = FALSE
) +
scale_fill_gradientn(colors = colors_cancer) +
guides(fill = guide_colorbar(title = "Activity")) +
ggtitle("ATAC Activity") # 添加标题
# RNA表达热图
DefaultAssay(sc) <- 'RNA'
p2 = DoHeatmap(sc,
features = top10$feature,
slot = "data",
draw.lines = TRUE,
raster = TRUE,
group.colors = identity_colors,
label = FALSE
) +
scale_fill_gradientn(colors = colors_cancer) +
guides(fill = guide_colorbar(title = "Expression")) +
ggtitle("RNA Expression") # 添加标题
# 组合并显示图形
combined_plot <- p1 + p2 # 使用/进行垂直排列
print(combined_plot)
#dir.create("_build/html/result/cluster_Diff", recursive = TRUE, showWarnings = FALSE)
# 保存图片
invisible(pdf("../result/cluster_Diff/Top10.GeneActivityandExpression.heatmap.pdf", width = 8, height = 10))
invisible(print(combined_plot))
invisible(dev.off())
}else{
print("没有基因活力显著差异的基因，请降低显著性阈值")
}
})
})小提琴图可视化差异结果
特征基因双模态小提琴图（ATAC_activity vs RNA）
可视化目的
- 对比各细胞亚群中高活性特征基因在染色质可及性（ATAC_activity）与转录表达（RNA）两种模态下的分布与一致性，辅助注释与机制解读。
图像构成
- 上半部分：ATAC 活性小提琴图，表示基因启动子/基因体的可及性分布（活力分数）。
- 下半部分：RNA 表达小提琴图，表示相同基因的转录水平分布。
- 选择与排序：每个细胞亚群选取基因活性最强的前 10 个特征基因；两种模态使用相同的基因集合与顺序；细胞按亚群分组并以颜色区分。
- 读图要点：小提琴的形状反映分布密度与异质性，颜色表示亚群，不同模态可直观比较同一基因在不同亚群的活性/表达强弱与分布差异。
解读要点
- 亚群特异性：目标亚群应呈现相应标志基因的活性/表达峰值与密度集中带。
- 跨模态一致性：ATAC 高活性与 RNA 高表达同步上升，支持开放染色质驱动的活跃转录。
- 不一致线索：活性高但表达低，可能提示转录预激活、时间滞后或远端增强子主导；表达高但活性低，可能受技术噪声或区域定义影响。
- 异质性识别：同一亚群内出现多峰/长尾，提示亚群结构或连续谱状态。
注意事项
- 间接性与噪声：基因活力为可及性的间接读数，较 RNA 更稀疏嘈杂，二者并非总一一对应。
R
# 先设置图形大小
suppressWarnings({
suppressMessages({
if(exists("top10")){
assays=c("ATAC_activity","RNA")
p=list()
for(assay in assays){
DefaultAssay(sc)=assay
if(length(top10$feature) > 1){
options(repr.plot.width = width, repr.plot.height = height)
stack=TRUE
p[[as.character(assay)]] = VlnPlot(sc,
features = top10$feature,
stack = stack,
pt.size = 0,
fill.by = "ident",
cols = identity_colors # 添加颜色参数
) +
xlab(if(assay == "ATAC_activity"){
"Gene Activity"
}else{
"Gene Expression"
}) +
theme_minimal() +
theme(
# 移除上方文本
axis.text.x.top = element_blank(),
strip.text.x = element_text(angle = 45, hjust = 0), # 分面标签角度# 设置下方基因名称
axis.text.x.bottom = element_text(
angle = 45,
hjust = 1, # 调整为1使文本右对齐
vjust = 1, # 调整为1使文本对齐轴线
face = "plain",
size = 8, # 添加字体大小设置
margin = margin(t = 0, r = 0, b = 0, l = 0)),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.ticks.x = element_blank(), # 移除刻度线
axis.line.x = element_blank(), # 移除轴线
plot.margin = margin(t = 5, r = 5, b = 5, l = 5) # 调整边距
) +NoLegend()+scale_x_discrete(labels = top10$feature)
}else if(length(top10$feature) ==1){
options(repr.plot.width = width, repr.plot.height = height)
stack=FALSE
p[[as.character(assay)]] = VlnPlot(sc,
features = top10$feature,
stack = stack,
pt.size = 0,
fill.by = "ident",
cols = identity_colors # 添加颜色参数
) +
xlab(if(assay == "ATAC_activity"){"Gene Activity"
}else{"Gene Expression"}) +
theme_minimal()
}
}
print(p[[1]])
print(p[[2]])
invisible(pdf("../result/cluster_Diff/Top10.GeneActivityandExpression.vlnplot.pdf",width=24,height=6))
invisible(print(p[[1]]))
invisible(print(p[[2]]))
invisible(dev.off())
}else{
print("没有基因活力显著差异的基因，请降低显著性阈值")
}
})
})气泡图可视化差异结果
双模态气泡图（ATAC_activity vs RNA）
可视化目的
- 并行展示各细胞类型中标志基因在染色质可及性与转录表达两种模态下的总体水平与检出情况，便于快速比较一致性与特异性。
图像构成与编码
- 维度：上方为 ATAC 活性气泡图（ATAC_activity），下方为 RNA 表达气泡图（RNA）。
- 行/列：行表示基因，列表示细胞类型/cluster（两张图使用相同的行列顺序以便对齐比较）。
- 气泡大小：该基因在对应细胞类型中被检测到的细胞比例（检测率；如活力/表达>0 的比例）。
- 颜色深浅：该基因在对应细胞类型的平均活力/平均表达水平（由浅到深表示由低到高）。
- 标准化建议：对两种模态分别在基因内做归一化/缩放，并统一配色范围与刻度标注，避免误读。
解读要点
- 特异性：目标细胞类型中应出现“大且深色”的气泡（高检测率且高均值），体现标志基因的特异活跃/高表达。
- 跨模态一致性：ATAC 与 RNA 在相同细胞类型同时出现“大且深色”的气泡，提示开放染色质与活跃转录的一致。
- 不一致线索：ATAC 大深而 RNA 小浅，可能为预激活/时间滞后或远端增强子主导；反之则可能受技术噪声或区域定义影响。
- 稀疏与异质：小但深的气泡可能来源于少数高值细胞；大但浅的气泡提示普遍低水平活跃/表达。
R
# 4. 气泡图
suppressWarnings({
suppressMessages({
if(exists("top10")){
options(repr.plot.width = width, repr.plot.height = height)
assays=c("ATAC_activity","RNA")
p=list()
for(assay in assays){
DefaultAssay(sc)=assay
p[[as.character(assay)]]=DotPlot(sc,
features = unique(top10$feature),
cols=c("#E5E3E3","#F0484A") # 从浅灰到深红的渐变
) +
#scale_size(range = c(0.5, 6)) + # 调整点的大小范围
# 修改颜色图例标题
guides(color = guide_colorbar(if(assay == "ATAC_activity"){title = "Average Activity"}else{title = "Average Expression"})) +
# 或者使用这种方式也可以
# labs(color = "Average Activity") +
theme(
axis.text.x.bottom = element_text(
angle = 45,
hjust = 1,
vjust = 1),
panel.grid.major = element_blank(), # 移除网格线
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"), # 添加坐标轴线
legend.position = "right" # 保留图例以查看效果
)
}
print(p[[1]])
print(p[[2]])
invisible(pdf("../result/cluster_Diff/Top10.GeneActivityandExpression.dotplot.pdf", width=24, height=6))
invisible(print(p[[1]]))
invisible(print(p[[2]]))
invisible(dev.off())
}else{
print("没有基因活力显著差异的基因，请降低显著性阈值")
}
})
})通过比较 scATAC-seq 数据计算的基因活力与 scRNA-seq 数据测得的基因表达水平，我们可以观察两者的一致性，这有助于我们理解染色质可及性与基因表达之间的潜在调控关系。同时，我们也可以通过观察经典 marker 基因的活性来验证细胞聚类的合理性。
但是，这些基因活力值的稳定性不如 scRNA-seq 直接测量的表达结果，主要有两个原因：
染色质数据本身具有较高的稀疏性；
基因调控区域的开放程度与基因的实际表达水平并非简单的线性关系。
尽管如此，这些基因活力数据仍然能够帮助我们较好地识别主要的细胞类型。但如果要将细胞进一步划分为更精细的亚群，仅依靠这种监督分析方法则较为困难。细胞分群间的相关性分析
跨细胞类型相关性分析（RNA 与 ATAC 活力）
分析目的
- 量化不同细胞类型在转录表达与染色质可及性两层面的相似性，检验聚类/注释合理性，并比较跨模态一致性与差异。
具体实现（与本代码一致）
- 输入：
sc（Seurat 对象）、cell_types_all（细胞类型列表）、celltype_colname（类型列名）。 - 数据来源：
sc@assays$RNA@data与sc@assays$ATAC_activity@data（均为对数标准化数据）。 - 类型汇总（pseudobulk）：
- 对每个细胞类型，取细胞层面的行均值：
Matrix::rowMeans(...)，生成基因×类型矩阵。 - 得到
rna_avg（尺寸约 25055 × n_types）与atac_avg（约 19620 × n_types）。
- 对每个细胞类型，取细胞层面的行均值：
- 相关性计算：
- 分别对
rna_avg与atac_avg按类型维度计算 Spearman 相关：cor(..., method = "spearman")。 - 输出两个类型×类型相关性矩阵：
RNA与ATAC，行列名均为cell_types_all。
- 分别对
- 返回结果：
list(RNA = rna_mat, ATAC = atac_mat)，用于热图可视化（附层次聚类与类型颜色条）。
- 输入：
解读要点
- 一致性：同谱系或功能相近的类型应在两张相关性热图中均表现出高相关。
- 跨模态对齐：同一类型在 RNA 与 ATAC 两个矩阵的相邻类型关系相似，支持注释可靠。
- 差异线索：RNA 高相关但 ATAC 低相关提示表观层面差异更明显或 ATAC 稀疏；反之可能反映预激活/时序差异。
- 异质性：相关性梯度或分支结构可提示亚群或连续谱状态。
R
# 计算不同细胞类型间的RNA和ATAC相关性
calculate_modality_correlations <- function(sc,cell_types,celltype_colname) {
# 初始化存储矩阵
n_types <- length(cell_types)
rna_mat <- matrix(0, n_types, n_types)
atac_mat <- matrix(0, n_types, n_types)
# 获取正确的基因数量
n_rna_genes <- nrow(sc@assays$RNA@data) # 应该是25055
n_atac_genes <- nrow(sc@assays$ATAC_activity@data) # 应该是19620
# 初始化正确维度的平均值矩阵
rna_avg <- matrix(0, n_rna_genes, n_types)
atac_avg <- matrix(0, n_atac_genes, n_types)
# 计算每个细胞类型的平均表达/可及性
for (idx in seq_along(cell_types)) {
current_type <- cell_types[idx]
# 获取当前类型的细胞索引
cell_idx <- which(sc@meta.data[[celltype_colname]] == current_type)
# 计算平均值
rna_avg[,idx] <- Matrix::rowMeans(sc@assays$RNA@data[, cell_idx])
atac_avg[,idx] <- Matrix::rowMeans(sc@assays$ATAC_activity@data[, cell_idx])
}
# 计算相关性
rna_mat <- cor(rna_avg, method = "spearman")
atac_mat <- cor(atac_avg, method = "spearman")
rownames(rna_mat) <- colnames(rna_mat) <- cell_types
rownames(atac_mat) <- colnames(atac_mat) <- cell_types
# 设置行列名
rownames(rna_mat) <- colnames(rna_mat) <- cell_types
rownames(atac_mat) <- colnames(atac_mat) <- cell_types
# 返回相关性矩阵
return(list(RNA = rna_mat, ATAC = atac_mat))
}
# 使用函数
results <- calculate_modality_correlations(sc,cell_types_all,celltype_colname)R
options(repr.plot.width = 18, repr.plot.height = 4)
# 定义通用可视化参数
get_breaks <- function(mat) {
min_break <- floor(min(mat) * 10) / 10 # 最小值向下取整到0.1
seq(min_break, 1, length.out = 101) # 固定最大值1
}
colors <- colorRampPalette(c("#4B2991", "#56B1F7", "#FCFFA4", "#F97B57", "#AB3282"))(100)
breaks <- get_breaks(c(results$RNA, results$ATAC)) # 统一计算色阶
# 通用热图参数
heatmap_config <- list(
color = colors,
breaks = breaks,
cluster_rows = FALSE,
cluster_cols = FALSE,
fontsize = 8,
border_color = NA,
silent = TRUE
)
# 生成热图
suppressWarnings({
rna_plot <- do.call(pheatmap, c(list(mat = results$RNA, main = "RNA Expression Correlation"), heatmap_config))
atac_plot <- do.call(pheatmap, c(list(mat = results$ATAC, main = "ATAC Accessibility Correlation"), heatmap_config))
# 展示并保存
grid.arrange(rna_plot$gtable, atac_plot$gtable, ncol = 2)
pdf("../result/modality_correlations.pdf", width = 12, height = 6)
grid.arrange(rna_plot$gtable, atac_plot$gtable, ncol = 2)
dev.off()
})* **左图：** RNA 表达相关性热图，反映了不同 cluster 在转录组水平上的相似程度和异质性
* **右图：** ATAC 可及性相关性热图，展示了染色质开放状态的相似性分布和异质性R
saveRDS(sc, file = "../result/result.rds")结果文件
├── cluster_Diff 细胞分群间基因活力差异分析结果文件夹
│ ├── cluster_Diff_genes_activity.txt 所有细胞分群间基因活力差异分析结果汇总表格
│ ├── FeaturePlot 降维图可视化每个细胞分群中差异倍数(logFC)top5 的基因活力情况
│ │ ├── feature_plot_cluster_ATAC_activity_.pdf
│ │ ├── feature_plot_cluster_ATAC_activity_.png
│ │ ├── feature_plot_cluster_RNA_.pdf
│ │ └── feature_plot_cluster_RNA_.png
│ ├── Top10.GeneActivityandExpression.dotplot.pdf 气泡图可视化每个细胞分群中差异倍数(logFC)top10 的基因活力和基因表达情况
│ ├── Top10.GeneActivityandExpression.heatmap.pdf 热图可视化每个细胞分群中差异倍数(logFC)top10 的基因活力和基因表达情况
│ └── Top10.GeneActivityandExpression.vlnplot.pdf 小提琴图可视化每个细胞分群中差异倍数(logFC)top10 的基因活力和基因表达情况
├── group_Diff 组间基因活力差异分析结果文件夹
│ ├── *.vs.heatmap.Activity.pdf 热图展示组间差异基因活力
│ ├── *.vs.heatmap.Activity.png 热图展示组间差异基因活力
│ ├── *.vs.heatmap.Expression.pdf 热图上述基因对应的基因表达情况
│ ├── *.vs.heatmap.Expression.png 热图上述基因对应的基因表达情况
│ ├── *.vs.Volcano.pdf 火山图可视化组间基因活力差异情况
│ ├── *.vs.Volcano.png 火山图可视化组间基因活力差异情况
│ └── group_Diff_genes_activity.txt 所有细胞分群中组间差异分析结果的汇总
└── modality_correlations.pdf
结果目录：目录下包括此项分析所涉及的所有图片以及表格等结果文件
文献案例解析
《Single-cell systems pharmacology identifies development-driven drug response and combination therapy in B cell acute lymphoblastic leukemia》
《Single-cell analyses reveal key immune cell subsets associated with response to PD-L1 blockade in triple-negative breast cancer》基于基因活力矩阵，对 scATAC-seq 数据进行细胞注释
《A cis-regulatory atlas in maize at single-cell resolution》
在这项研究中，研究人员通过整合单细胞水平的 ATAC-seq 和 RNA-seq 数据，探究了玉米基因组中基因可及性与基因表达之间的关系。研究发现，在 36,322 个分析的基因中，大多数基因表现出基因可及性与转录活性的正相关关系（图 E、F、G、H）。然而，研究也发现约 6,063 个基因虽然具有较高的染色质可及性，但却没有显著的转录活性（图 H），这表明基因的可及性并不必然导致基因表达。为了解释这一现象，研究人员对基因进行了表观遗传学分析。结果表明，这些可及但不表达的基因富集了 H3K27me3 修饰（图 J）；相比之下，那些既不可及又不表达的基因（约 4,315 个）则主要通过 DNA 甲基化来维持其沉默状态。这一发现揭示了玉米基因沉默调控的两种主要机制，（1）通过 H3K27me3 修饰，对基因进行沉默；（2）依赖 DNA 甲基化来维持沉默状态
参考资料
[1] Stuart T, Srivastava A, Madad S, et al.Single-cell chromatin state analysis with Signac (vol 18, pg 1333, 2021)[J].Nature methods, 2022(2):19. [2] Hao Y, Stuart T, Kowalski M, et al.Dictionary learning for integrative, multimodal, and scalable single-cell analysis[J].bioRxiv, 2022.DOI: 10.1101/2022.02.24.481684.
附录
.txt ：结果数据表格文件，文件以制表符（Tab）分隔。unix/Linux/Mac 用户使用 less 或 more 命令查看；windows 用户使用高级文本编辑器 Notepad++ 等查看，也可以用 Microsoft Excel 打开。
.pdf ：结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用 Adobe Illustrator 进行图片编辑，用于文章发表等。
.rds : 包含基因活力 assay 的 Seurat 对象，需要在 R 环境中打开查看和用于进一步分析。