在本教程中，我们将学习使用Signac包进行DNA序列的motif富集分析。Signac包可以使用两种互补的方法进行motif分析：一种是通过在一组差异可及性peaks中找到overrepresented的基序motifs，另一种是在不同细胞组之间执行差异基序活性（motif activity）分析的方法。

安装并加载所需的R包

  1. if (!requireNamespace("BiocManager", quietly = TRUE))
  2.     install.packages("BiocManager")
  3. BiocManager::install("JASPAR2018")
  4. BiocManager::install("TFBSTools")
  6. library(Signac)
  7. library(Seurat)
  8. library(JASPAR2018)
  9. library(TFBSTools)
  10. library(BSgenome.Mmusculus.UCSC.mm10)
  11. library(patchwork)
  12. set.seed(1234)

加载所需的数据集

  1. mouse_brain <- readRDS("../vignette_data/adult_mouse_brain.rds")
  2. mouse_brain
  3. ## An object of class Seurat 
  4. ## 178653 features across 3503 samples within 2 assays 
  5. ## Active assay: peaks (157203 features, 157203 variable features)
  6. ##  1 other assay present: RNA
  7. ##  2 dimensional reductions calculated: lsi, umap
  9. p1 <- DimPlot(mouse_brain, label = TRUE, pt.size = 0.1) + NoLegend()
  10. p1

使用Signac包进行单细胞ATAC-seq数据分析（二）：Motif analysis with Signac - 图1

构建Motif类

为了有助于Signac进行motif富集分析，我们需要创建一个Motif类来存储所有必需的信息，其中包括位置权重矩阵（PWM）或位置频率矩阵（PFM）的列表以及motif出现矩阵。在这里，我们构造了一个Motif对象，并将其添加到我们的小鼠大脑数据集中。可以使用AddMotifObject函数将motif对象添加到Seurat对象的assay中：

  1. # Get a list of motif position frequency matrices from the JASPAR database
  2. # 使用getMatrixSet函数从JASPAR数据库中提取Motif的PFM矩阵信息
  3. pfm <- getMatrixSet(
  4.   x = JASPAR2018,
  5.   opts = list(species = 9606, all_versions = FALSE)
  6. )
  8. # Scan the DNA sequence of each peak for the presence of each motif
  9. # 使用CreateMotifMatrix函数构建Motif矩阵对象
  10. motif.matrix <- CreateMotifMatrix(
  11.   features = StringToGRanges(rownames(mouse_brain), sep = c(":", "-")),
  12.   pwm = pfm,
  13.   genome = 'mm10',
  14.   sep = c(":", "-"),
  15.   use.counts = FALSE
  16. )
  18. # Create a new Mofif object to store the results
  19. # 使用CreateMotifObject函数构建Motif对象
  20. motif <- CreateMotifObject(
  21.   data = motif.matrix,
  22.   pwm = pfm
  23. )
  25. # Add the Motif object to the assay
  26. # 使用AddMotifObject函数将Motif类添加到Seurat对象中
  27. mouse_brain[['peaks']] <- AddMotifObject(
  28.   object = mouse_brain[['peaks']],
  29.   motif.object = motif
  30. )

为了找出overrepresented的motifs基序，我们还需要计算peaks的一些序列特征，例如GC含量，序列长度和二核苷酸频率。我们可以使用RegionStats函数计算这些信息，并将结果存储在Seurat对象的元数据中。

  1. # 使用RegionStats函数计算peaks区域的序列特生
  2. mouse_brain <- RegionStats(
  3.   object = mouse_brain,
  4.   genome = BSgenome.Mmusculus.UCSC.mm10,
  5.   sep = c(":", "-")
  6. )

Finding overrepresented motifs

为了鉴定出潜在的重要细胞类型特异性调控序列，我们可以富集出在不同细胞类型之间差异可及性的peaks中overrepresented的DNA基序。
这里，我们首先鉴定出Pvalb和Sst细胞类群之间的差异可及性peaks。然后，对它们进行一次超几何测试检验，以检验在给定频率下偶然观察到基序的可能性，并与GC含量匹配的背景峰进行比较。

  1. # 使用FindMarkers函数鉴定差异可及性peaks
  2. da_peaks <- FindMarkers(
  3.   object = mouse_brain,
  4.   ident.1 = 'Pvalb',
  5.   ident.2 = 'Sst',
  6.   only.pos = TRUE,
  7.   test.use = 'LR',
  8.   latent.vars = 'nCount_peaks'
  9. )
  11. # Test the differentially accessible peaks for overrepresented motifs
  12. # 使用FindMotifs函数进行motif富集分析
  13. enriched.motifs <- FindMotifs(
  14.   object = mouse_brain,
  15.   features = head(rownames(da_peaks), 1000)
  16. )
  18. # 查看motif富集分析的结果
  19. # sort by p-value and fold change
  20. enriched.motifs <- enriched.motifs[order(enriched.motifs[, 7], -enriched.motifs[, 6]), ]
  21. head(enriched.motifs)
  22. ##             motif observed background percent.observed percent.background
  23. ## MA0497.1 MA0497.1      431       9019             43.1            22.5475
  24. ## MA1151.1 MA1151.1      226       3227             22.6             8.0675
  25. ## MA0773.1 MA0773.1      304       5421             30.4            13.5525
  26. ## MA0072.1 MA0072.1      218       3160             21.8             7.9000
  27. ## MA0052.3 MA0052.3      315       5844             31.5            14.6100
  28. ## MA1150.1 MA1150.1      211       3108             21.1             7.7700
  29. ##          fold.enrichment       pvalue  motif.name
  30. ## MA0497.1        1.911520 1.549646e-48       MEF2C
  31. ## MA1151.1        2.801363 5.818264e-47        RORC
  32. ## MA0773.1        2.243129 1.365161e-44       MEF2D
  33. ## MA0072.1        2.759494 4.091607e-44 RORA(var.2)
  34. ## MA0052.3        2.156057 6.539737e-43       MEF2A
  35. ## MA1150.1        2.715573 1.600301e-41        RORB

我们还可以使用MotifPlot函数绘制motif的位置权重矩阵，可视化不同的motif序列。

  1. MotifPlot(
  2.   object = mouse_brain,
  3.   motifs = head(rownames(enriched.motifs))
  4. )

使用Signac包进行单细胞ATAC-seq数据分析（二）：Motif analysis with Signac - 图2

Computing motif activities

我们还可以通过运行chromVAR计算每个细胞的基序活性得分(motif activity score)，这样我们可以查看每个细胞的motif activities，并且还提供了一种识别细胞类型之间差异激活的基序的替代方法。ChromVAR可识别与细胞之间染色质可及性变异相关的基序，有关该方法的完整说明，请参见chromVAR的说明文档。

  1. # 使用RunChromVAR函数计算所有细胞中的motif activities
  2. mouse_brain <- RunChromVAR(
  3.   object = mouse_brain,
  4.   genome = BSgenome.Mmusculus.UCSC.mm10
  5. )
  7. DefaultAssay(mouse_brain) <- 'chromvar'
  9. # look at the activity of Mef2c, the top hit from the overrepresentation testing
  10. p2 <- FeaturePlot(
  11.   object = mouse_brain,
  12.   features = rownames(enriched.motifs)[[1]],
  13.   min.cutoff = 'q10',
  14.   max.cutoff = 'q90',
  15.   pt.size = 0.1
  16. )
  17. p1 + p2

使用Signac包进行单细胞ATAC-seq数据分析（二）：Motif analysis with Signac - 图3

我们还可以直接测试不同细胞类型之间motif的差异活性得分，这和对不同细胞类型之间的差异可及性peaks进行富集测试的结果相类似。

  1. differential.activity <- FindMarkers(
  2.   object = mouse_brain,
  3.   ident.1 = 'Pvalb',
  4.   ident.2 = 'Sst',
  5.   only.pos = TRUE,
  6.   test.use = 'LR',
  7.   latent.vars = 'nCount_peaks'
  8. )
  10. # 使用MotifPlot函数对富集到的motif进行可视化
  11. MotifPlot(
  12.   object = mouse_brain,
  13.   motifs = head(rownames(differential.activity)),
  14.   assay = 'peaks'
  15. )

使用Signac包进行单细胞ATAC-seq数据分析（二）：Motif analysis with Signac - 图4

参考来源：https://satijalab.org/signac/articles/motif_vignette.html