建议开始此流程前先学习:
注意:该流程某些步骤需要较优的CPU和较大的内存,建议在集群上运行
亲测,CPU为7840h,内存为40G(所谓的运行内存)会在第四步运行失败,在集群中使用80G内存则运行成功
在集群中运行时请善用 ggsave() 和 saveRDS() 函数保存所需图片和数据
一、安装所需的包和数据
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("JASPAR2020") # JASPAR2020 数据包包含了 JASPAR 数据库 2020 版的转录因子结合位点信息,包括各种生物物种中的转录因子结合位点的位置频率矩阵、Motif 描述、注释等信息
BiocManager::install("TFBSTools")
BiocManager::install("BSgenome.Mmusculus.UCSC.mm10") # 这个数据包包含了小鼠基因组(mm10 版本)的信息
BiocManager::install("motifmatchr")
BiocManager::install("chromVAR")
install.packages('ggseqlogo')
install.packages("patchwork")
install.packages("Signac")
install.packages("Seurat")
install.packages("ggplot2")
library(Signac)
library(Seurat)
library(JASPAR2020)
library(motifmatchr)
library(ggseqlogo)
library(chromVAR)
library(TFBSTools)
library(BSgenome.Mmusculus.UCSC.mm10)
library(patchwork)
library(ggplot2)
set.seed(1234) # 这一行代码设置随机数种子,以便结果可以重复复现
Analyzing adult mouse brain scATAC-seq • Signac (stuartlab.org)https://stuartlab.org/signac/articles/mouse_brain_vignette adult_mouse_brain.rds 这一数据集是由上方网址中的 scATAC-seq 分析流程得出(具体步骤可参考: scATAC-seq 数据分析: Signac(官方流程复现 + 步骤详解)_蚁丘的博客-CSDN博客,注意,这与官方使用的数据不一样,只供参考其流程),在最后一步完成后,运行下方代码即可得到
saveRDS(object = brain, file = "../vignette_data/adult_mouse_brain.rds")
方便起见,本文直接给出 adult_mouse_brain.rds 的下载地址,地址如下:https://pan.baidu.com/s/1WuNLK5PY2NcjvssIBzgzMw?pwd=o0xe https://pan.baidu.com/s/1WuNLK5PY2NcjvssIBzgzMw?pwd=o0xe%C2%A0 请注意!由于本文作者在得到 adult_mouse_brain.rds 文件的流程中,使用的原始数据并未更新,如果直接使用给出的文件进行分析,得出的结果会与官方最新参考流程的结果有所偏差,本文所使用的为本文作者提供的数据。
# 加载所需数据集
mouse_brain <- readRDS("../vignette_data/adult_mouse_brain.rds")
mouse_brain
# An object of class Seurat
# 298331 features across 3517 samples within 2 assays
# Active assay: peaks (276523 features, 276523 variable features)
# 1 other assay present: RNA
# 2 dimensional reductions calculated: lsi, umap
p1 <- DimPlot(mouse_brain, label = TRUE, pt.size = 0.1) + NoLegend()
p1
二、为 Seurat 对象添加 motif 信息
通过以下代码,可从 JASPAR 数据库获取的转录因子结合位点信息与单细胞数据集进行关联,从而为之后的motif分析提供必要的数据基础。这有助于将转录因子调控信息与单细胞染色质分析相结合,深入研究细胞的转录调控网络。
# 从JASPAR数据库获取转录因子结合位点的位置频率矩阵
pfm <- getMatrixSet(
x = JASPAR2020, # 使用JASPAR2020版本的数据库
opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)
# collection = "CORE":指定要获取的Motif集合。在这里,选择了JASPAR数据库中的"CORE"集合,该集合包含了基础、常见的转录因子结合位点Motif信息。
# tax_group = 'vertebrates':指定生物分类,以便获取特定生物类群的Motif信息。在这里,选择了"vertebrates"(脊椎动物)
# 将获取的Motif信息添加到小鼠脑数据集中
mouse_brain <- AddMotifs(
object = mouse_brain, # 要添加到的Seurat对象
genome = BSgenome.Mmusculus.UCSC.mm10, # 对应的小鼠基因组信息
pfm = pfm # 之前获取的位置频率矩阵
)
三、寻找过表达的 motif
-
鉴定细胞类型特异性调控序列: 通过寻找在不同细胞类型之间可及性存在差异的峰集中的DNA motifs,可以鉴定可能的细胞类型特异性调控序列。这些 motifs 可能代表着在某些细胞类型中特定的转录因子结合位点。
-
在细胞类型之间找到差异可及性峰集: 在 Pvalb 和 Sst 抑制性间脑内神经元之间寻找差异可及性峰集。既寻找在这两种细胞类型中表现出显著差异的开放染色质区域。
-
调整 FindMarkers() 函数中的 min.pct 阈值: 对于稀疏数据(如 scATAC-seq 数据),发现在 FindMarkers() 函数中降低 min.pct 阈值通常是必要的。min.pct 阈值用于控制用于寻找差异峰的最小细胞百分比阈值。作者指出,对于 scATAC-seq 数据,需要将默认值 0.1 降低,因为这个阈值最初是为了处理 scRNA-seq 数据设计的。
-
执行超几何检验: 为了测试在给定频率下观察到的 motif 是否可能是偶然发生的,需要执行了超几何检验。超几何检验用于检查 motif 在给定频率下是否与背景峰集(具有匹配 GC 含量的峰集)中的出现概率相比呈现显著差异。
寻找差异可及性峰:
# 使用 FindMarkers() 函数寻找差异可及性峰集
da_peaks <- FindMarkers(
object = mouse_brain, # 使用的 Seurat 数据对象
ident.1 = 'Pvalb', # 第一个细胞群的标识
ident.2 = 'Sst', # 第二个细胞群的标识
only.pos = TRUE, # 只考虑正的 log fold change
test.use = 'LR', # 使用 Likelihood Ratio test 进行差异分析
min.pct = 0.05, # 最小细胞百分比阈值(适用于 scATAC-seq 数据)
latent.vars = 'nCount_peaks' # 指定隐变量,这里使用峰的计数
)
# 从差异分析结果中获取显著差异的可及性峰集
top.da.peak <- rownames(da_peaks[da_peaks$p_val < 0.005, ])
下方代码使用 FindMotifs()
函数来测试上一步找到的差异可及性峰集是否富集了特定的 DNA motifs。函数将分析这些峰集中是否有与已知转录因子结合位点相关的 DNA motifs 出现的富集情况:
# 使用 FindMotifs() 函数进行富集分析
enriched.motifs <- FindMotifs(
object = mouse_brain,
features = top.da.peak
)
# Selecting background regions to match input sequence characteristics
# Matching GC.percent distribution
# Testing motif enrichment in 3805 regions
View(head(enriched.motifs))
motif
: DNA motif 的标识observed
: 在差异可及性峰集中观察到的该 DNA motif 的次数background
: 在背景峰集中出现的该 DNA motif 的次数percent.observed
: 该 DNA motif 在差异可及性峰集中的百分比percent.background
: 该 DNA motif 在背景峰集中的百分比fold.enrichment
: 富集倍数,即观察到的百分比除以背景的百分比pvalue
: p 值,表示观察到的富集是否显著motif.name
: DNA motif 的名称p.adjust
: 经过多重检验校正的调整 p 值
接下来使用 MotifPlot()对富集的 DNA motifs 进行可视化:
MotifPlot(
object = mouse_brain,
motifs = head(rownames(enriched.motifs))
)
四、计算 motif 的活性
通过运行 chromVAR()来计算每个细胞的 motif 活性分数。计算出分数后,我们能够可视化每个细胞的 motif 活性
请注意!!!RunChromVAR() 函数运行所占CPU和内存容量大,推荐在集群中运行下述代码(亲测,CPU为7840h,内存为40G会运行失败,在集群中使用80G内存则运行成功)
# 该步骤占用内存和CPU大,推荐在集群中运行
# 运行 ChromVAR,计算每个细胞的 Motif 活性分数
mouse_brain <- RunChromVAR(
object = mouse_brain,
genome = BSgenome.Mmusculus.UCSC.mm10 # 使用的基因组信息
)
# 将计算得到的 Motif 活性分数设置为 Seurat 对象的默认分析集
DefaultAssay(mouse_brain) <- 'chromvar'
# 创建一个特征图,展示 Motif "MA0497.1" 的活性情况
p2 <- FeaturePlot(
object = mouse_brain,
features = "MA0497.1", # 要可视化的 Motif 名称
min.cutoff = 'q10', # 活性值的最小可视化范围
max.cutoff = 'q90', # 活性值的最大可视化范围
pt.size = 0.1
)
p1 + p2
下段代码通过 FindMarkers()
函数计算两个细胞类型之间的差异活性分数,并使用 MotifPlot()
函数绘制与这些差异活性分数相关的基因组特征的 motif 图,从而帮助分析两个细胞类型之间的转录因子活性差异
# 使用FindMarkers函数计算两个细胞类型之间的差异活性分数
differential.activity <- FindMarkers(
object = mouse_brain,
ident.1 = 'Pvalb', # 第一个细胞类型的标识符
ident.2 = 'Sst', # 第二个细胞类型的标识符
only.pos = TRUE, # 仅考虑正向的活性差异
mean.fxn = rowMeans, # 在计算折叠变化时使用行平均值函数
fc.name = "avg_diff" # 折叠变化值的名称
)
# 使用MotifPlot函数绘制与差异活性分数相关的基因组特征的模体图
MotifPlot(
object = mouse_brain,
motifs = head(rownames(differential.activity)),
assay = 'peaks' # 使用基因组特征的'peaks'数据
)