使用 Signac 进行 TF motif 分析

蚁丘

已于 2023-08-28 16:45:13 修改

阅读量511

点赞数 1

文章标签： r语言数据分析

于 2023-08-28 16:41:56 首次发布

本文链接：https://blog.csdn.net/ZIGRA/article/details/132362313

版权

建议开始此流程前先学习：

scATAC-seq 数据分析: Signac（官方流程复现 + 步骤详解）_蚁丘的博客-CSDN博客https://blog.csdn.net/ZIGRA/article/details/132312220?spm=1001.2014.3001.5501

注意：该流程某些步骤需要较优的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'数据
)

参考文章：
Motif analysis with Signac • Signac (stuartlab.org)https://stuartlab.org/signac/articles/motif_vignette

蚁丘

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
使用 Signac 进行 TF motif 分析

scATAC-seq 数据分析: Signac（官方流程复现 + 步骤详解）_蚁丘的博客-CSDN博客亲测，CPU为7840h，内存为40G（所谓的运行内存）会在第四步运行失败，在集群中使用80G内存则运行成功在集群中运行时请善用 ggsave() 和 saveRDS() 函数保存所需图片和数据。
复制链接

扫一扫