单细胞 scATAC-seq & scRNA-seq 双组学标准分析：GeneActivity 分析
文档概述
GeneActivity 分析通过计算基因体及其上游调控区域的染色质可及性，推断潜在的转录活性，以弥补单细胞 ATAC 数据中“只有开放区域，没有基因表达量”的信息缺口。本技术文档围绕核心概念、计算方法、可视化解读、云平台操作及常见问题展开，帮助您快速理解基因活力分析的完整过程及相关结果。
一、核心概念
1. 什么是 GeneActivity？
- 基因活力（Gene Activity）：统计每个基因的启动子（通常取 TSS 上游 2 kb）与基因体区域内的 fragment 数量，经过标准化后得到的活性值。该值反映“该基因是否具备转录潜力”，而非实际 mRNA 表达量。
- 片段（fragment）：Tn5 转座酶在开放染色质区域切割所得的 DNA 片段，是 ATAC 信号的载体。
- 调控区域覆盖：默认采用“基因体 + 上游 2 kb”框定统计范围，既包含核心启动子也记录转录延伸阶段的可及性。
2. 基因活力 vs. 基因表达
| 指标 | 数据来源 | 反映信息 | 常见用途 |
|---|---|---|---|
| Gene Activity | scATAC-seq | 染色质开放程度 → 潜在转录能力 | 细胞类型注释、差异可及性基因、跨组学比对 |
| Gene Expression | scRNA-seq | mRNA 实际丰度 | 传统差异表达、功能富集、轨迹分析 |
TIP
GeneActivity 与 Gene Expression 通常正相关，但也会因表观遗传沉默、转录后调控或技术噪音而出现不一致。对比两者可更深入理解调控机制。
3. GeneActivity 值的表示
GeneActivity 输出为稀疏矩阵，行对应基因、列对应单细胞，矩阵元素为基因区域的 fragment 计数，例如：
6 x 3 sparse Matrix of class "dgCMatrix"
AAACCCAATGCGTGGCA_1 AAACCCAATGTGCTACT_1 AAACCCACGACATCCGG_1
Hnf4g . . .
Zfhx4 . 5 2
Pex2 . . 2
UBC . . .
1700008P02Rik . . .
Pkia . . .二、GeneActivity 的计算与差异分析
1. 计算方法
以 Signac 为例，GeneActivity 分析分为三个阶段：定义基因区域、统计 fragment、标准化输出。典型代码如下：
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
# 1. 构建染色质 assay 并添加注释
annotation <- GetGRangesFromEnsDb(EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"
genome(annotation) <- "hg38"
Annotation(atac_obj) <- annotation
# 2. 计算基因活力矩阵（假设 atac_obj 已经完成质控、标准化）
gene.activities <- GeneActivity(
object = atac_obj,
extend.upstream = 2000, # 上游 2kb
extend.downstream = 0
)
# 3. 创建 GeneActivity assay 并进入标准化流程
atac_obj[["GeneActivity"]] <- CreateAssayObject(counts = gene.activities)
DefaultAssay(atac_obj) <- "GeneActivity"
atac_obj <- NormalizeData(atac_obj)
atac_obj <- ScaleData(atac_obj)TIP
extend.upstream参数可根据物种或研究对象调整（1–3 kb 最常见）。- 计算 GeneActivity 之前需确保 fragment 文件与细胞条形码对齐，并完成基础 QC。
2. 差异 GeneActivity 计算
GeneActivity assay 构建完成后，可直接复用 FindMarkers() 或 FindAllMarkers() 完成差异分析，与 RNA 数据的差异基因分析方法一样：
DefaultAssay(atac_obj) <- "GeneActivity"
gene_activity_markers <- FindAllMarkers(
object = atac_obj,
group.by = "seurat_cluster",
min.pct = 0.1,
logfc.threshold = 0.25,
only.pos = TRUE,
test.use = "wilcox"
)- 适用场景：一次性比较所有细胞群体（
FindAllMarkers），或指定两组（FindMarkers）完成疾病/处理 vs 对照的对比。 - 推荐筛选：
p_val_adj < 0.05、avg_log2FC > 0.25、pct_in > pct_out。
3. 结果字段说明
GeneActivity 差异表格与差异 Gene 表类似，包含：
feature：基因名称；cluster/group：对应细胞类型或分组标签；avgExpr/logFC：平均活力及 log2FC；pct_in/pct_out：在目标组和背景组的检出率；statistic、auc、p_val、p_val_adj：统计显著性指标。
三、GeneActivity 结果展示
1. 差异表格
交互式表格功能强大，可根据 logFC、p 值及检出率等多维指标灵活筛选目标基因，支持一键导出结果，方便后续进行注释、功能富集分析或个性化可视化展示。
云平台默认根据 logFC Top基因生成多种可视化，串联 ATAC 与 RNA 两个模态，帮助快速定位核心调控事件。
2. FeaturePlot
利用FeaturePlot 函数同步展示 ATAC 基因活力与 RNA 表达，便于一图对比染色质开放和转录水平。图中每个点代表一个细胞，颜色深浅分别反映目标基因在 ATAC（活力）和 RNA（表达）两个层面的丰度。通过该图，可以直观验证核心 Top 差异基因是否在特定细胞群体中既开放又表达，实现多组学结果的快速交叉验证。
3. 双热图
左侧展示 GeneActivity（ATAC 活力）热图，右侧为 RNA 表达热图。每一行为核心 Top 差异基因，每一列对应单个细胞，顶部通过不同颜色表示细胞分群。通过两者对比，能快速评估染色质开放水平与真实转录表达的一致性及潜在差异，有助于识别仅开放未表达、仅表达未开放等特殊调控现象。
4. 小提琴图
每个基因对应一组 ATAC（基因活力）和 RNA（基因表达）的小提琴图，细致展示各细胞群体在不同模态下的分布特征。该图有助于快速定位高活力但低表达、或高表达但活力较低的基因，揭示潜在的表观遗传调控机制与转录后调控现象。
5. 气泡图
双层气泡图通过气泡大小直观反映检出率（pct_in），颜色深浅代表平均基因活力或表达量。此图可在细胞群体与基因的矩阵中高效定位显著差异基因，一目了然突出核心调控特征。
TIP：本节所有可视化（FeaturePlot降维图、双热图、小提琴图、气泡图）均为“基因活力差异倍数Top基因”在RNA层面（表达量）与ATAC层面（基因活力）的对比展示。通过双模态对照，能更直观地识别核心调控基因及其跨层面的表达开放关系。
四、多组学相关性分析
本部分的主要目的是比较不同细胞群（类型）在RNA表达和ATAC基因可及性两个层面的异质性与相关性。具体流程为：首先分别计算每个细胞群的RNA平均表达量和ATAC基因活力均值，生成以细胞群为单位的均值矩阵；然后基于这些均值矩阵，采用Spearman秩相关系数，对所有细胞群进行两两相关性分析，分别得到RNA表达相关性矩阵和ATAC活力相关性矩阵。
- RNA相关性矩阵：评估细胞群间转录表达模式的异质性与相似性；
- ATAC活力相关性矩阵：反映各细胞群在染色质可及性（开放状态）上的一致性或多样性。
相关性热图的左侧为各细胞群（或类型）间的RNA基因表达相关性，右侧为ATAC基因活力相关性。横轴和纵轴均为细胞群，每个方格的颜色代表两个细胞群体的相关系数（表达或活力相关性）：颜色越红相关性越高，相似性更强；颜色越蓝则相关性低、异质性大。通过热图对比，可直观评估哪些细胞群体在RNA和ATAC层面均高度相关，哪些则只在某一层面相关，从而揭示表观遗传调控和转录表达之间的一致性或解耦现象。
解读建议：
- 通过比较相关性热图，整体观察ATAC活力与RNA表达的相关性分布。如果ATAC层面的相关性普遍高于RNA，说明表达的异质性大于基因活力的异质性；反之，若RNA相关性更高，则基因活力异质性相对更大。
- 通常在同一套单细胞多组学数据中，虽然ATAC peaks的开放常暗示更强的异质性，但以GeneActivity为代表的ATAC层面活力的异质性通常小于RNA表达的异质性。具体需结合热图和数据实际情况，深入分析不同层面对细胞群异质性的生物学意义。
五、常见问题
Q1：GeneActivity 值为何与 RNA 表达不完全一致？
A：染色质开放只代表具备转录条件，仍可能受到抑制性组蛋白修饰、缺乏关键转录因子或转录后调控的影响。相反，RNA 高但 GeneActivity 低多由 ATAC 稀疏或远程增强子调控造成。
Q2：为何默认采用“基因体 + 上游 2 kb”？可否调整？
A：该范围兼顾启动子与转录延伸阶段，既具生物学意义又能缓解 ATAC 数据稀疏问题。extend.upstream 参数可按需调整（如 1 kb / 2.5 kb），但需保持全数据分析的一致性。