在本教程中，我們將學習使用Signac包進行DNA序列的motif富集分析。Signac包可以使用兩種互補的方法進行motif分析：一種是通過在一組差異可及性peaks中找到overrepresented的基序motifs，另一種是在不同細胞組之間執(zhí)行差異基序活性（motif activity）分析的方法。

安裝并加載所需的R包

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("JASPAR2018")
BiocManager::install("TFBSTools")

library(Signac)
library(Seurat)
library(JASPAR2018)
library(TFBSTools)
library(BSgenome.Mmusculus.UCSC.mm10)
library(patchwork)
set.seed(1234)

加載所需的數據集

mouse_brain <- readRDS("../vignette_data/adult_mouse_brain.rds")
mouse_brain
## An object of class Seurat 
## 178653 features across 3503 samples within 2 assays 
## Active assay: peaks (157203 features, 157203 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

image

構建Motif類

為了有助于Signac進行motif富集分析，我們需要創(chuàng)建一個Motif類來存儲所有必需的信息，其中包括位置權重矩陣（PWM）或位置頻率矩陣（PFM）的列表以及motif出現(xiàn)矩陣。在這里，我們構造了一個Motif對象，并將其添加到我們的小鼠大腦數據集中?？梢允褂?code>AddMotifObject函數將motif對象添加到Seurat對象的assay中：

# Get a list of motif position frequency matrices from the JASPAR database
# 使用getMatrixSet函數從JASPAR數據庫中提取Motif的PFM矩陣信息
pfm <- getMatrixSet(
  x = JASPAR2018,
  opts = list(species = 9606, all_versions = FALSE)
)

# Scan the DNA sequence of each peak for the presence of each motif
# 使用CreateMotifMatrix函數構建Motif矩陣對象
motif.matrix <- CreateMotifMatrix(
  features = StringToGRanges(rownames(mouse_brain), sep = c(":", "-")),
  pwm = pfm,
  genome = 'mm10',
  sep = c(":", "-"),
  use.counts = FALSE
)

# Create a new Mofif object to store the results
# 使用CreateMotifObject函數構建Motif對象
motif <- CreateMotifObject(
  data = motif.matrix,
  pwm = pfm
)

# Add the Motif object to the assay
# 使用AddMotifObject函數將Motif類添加到Seurat對象中
mouse_brain[['peaks']] <- AddMotifObject(
  object = mouse_brain[['peaks']],
  motif.object = motif
)

為了找出overrepresented的motifs基序，我們還需要計算peaks的一些序列特征，例如GC含量，序列長度和二核苷酸頻率。我們可以使用RegionStats函數計算這些信息，并將結果存儲在Seurat對象的元數據中。

# 使用RegionStats函數計算peaks區(qū)域的序列特生
mouse_brain <- RegionStats(
  object = mouse_brain,
  genome = BSgenome.Mmusculus.UCSC.mm10,
  sep = c(":", "-")
)

Finding overrepresented motifs

為了鑒定出潛在的重要細胞類型特異性調控序列，我們可以富集出在不同細胞類型之間差異可及性的peaks中overrepresented的DNA基序。
這里，我們首先鑒定出Pvalb和Sst細胞類群之間的差異可及性peaks。然后，對它們進行一次超幾何測試檢驗，以檢驗在給定頻率下偶然觀察到基序的可能性，并與GC含量匹配的背景峰進行比較。

# 使用FindMarkers函數鑒定差異可及性peaks
da_peaks <- FindMarkers(
  object = mouse_brain,
  ident.1 = 'Pvalb',
  ident.2 = 'Sst',
  only.pos = TRUE,
  test.use = 'LR',
  latent.vars = 'nCount_peaks'
)

# Test the differentially accessible peaks for overrepresented motifs
# 使用FindMotifs函數進行motif富集分析
enriched.motifs <- FindMotifs(
  object = mouse_brain,
  features = head(rownames(da_peaks), 1000)
)

# 查看motif富集分析的結果
# sort by p-value and fold change
enriched.motifs <- enriched.motifs[order(enriched.motifs[, 7], -enriched.motifs[, 6]), ]
head(enriched.motifs)
##             motif observed background percent.observed percent.background
## MA0497.1 MA0497.1      431       9019             43.1            22.5475
## MA1151.1 MA1151.1      226       3227             22.6             8.0675
## MA0773.1 MA0773.1      304       5421             30.4            13.5525
## MA0072.1 MA0072.1      218       3160             21.8             7.9000
## MA0052.3 MA0052.3      315       5844             31.5            14.6100
## MA1150.1 MA1150.1      211       3108             21.1             7.7700
##          fold.enrichment       pvalue  motif.name
## MA0497.1        1.911520 1.549646e-48       MEF2C
## MA1151.1        2.801363 5.818264e-47        RORC
## MA0773.1        2.243129 1.365161e-44       MEF2D
## MA0072.1        2.759494 4.091607e-44 RORA(var.2)
## MA0052.3        2.156057 6.539737e-43       MEF2A
## MA1150.1        2.715573 1.600301e-41        RORB

我們還可以使用MotifPlot函數繪制motif的位置權重矩陣，可視化不同的motif序列。

MotifPlot(
  object = mouse_brain,
  motifs = head(rownames(enriched.motifs))
)

image

Computing motif activities

我們還可以通過運行chromVAR計算每個細胞的基序活性得分(motif activity score)，這樣我們可以查看每個細胞的motif activities，并且還提供了一種識別細胞類型之間差異激活的基序的替代方法。ChromVAR可識別與細胞之間染色質可及性變異相關的基序，有關該方法的完整說明，請參見chromVAR的說明文檔。

# 使用RunChromVAR函數計算所有細胞中的motif activities
mouse_brain <- RunChromVAR(
  object = mouse_brain,
  genome = BSgenome.Mmusculus.UCSC.mm10
)

DefaultAssay(mouse_brain) <- 'chromvar'

# look at the activity of Mef2c, the top hit from the overrepresentation testing
p2 <- FeaturePlot(
  object = mouse_brain,
  features = rownames(enriched.motifs)[[1]],
  min.cutoff = 'q10',
  max.cutoff = 'q90',
  pt.size = 0.1
)
p1 + p2

image

我們還可以直接測試不同細胞類型之間motif的差異活性得分，這和對不同細胞類型之間的差異可及性peaks進行富集測試的結果相類似。

differential.activity <- FindMarkers(
  object = mouse_brain,
  ident.1 = 'Pvalb',
  ident.2 = 'Sst',
  only.pos = TRUE,
  test.use = 'LR',
  latent.vars = 'nCount_peaks'
)

# 使用MotifPlot函數對富集到的motif進行可視化
MotifPlot(
  object = mouse_brain,
  motifs = head(rownames(differential.activity)),
  assay = 'peaks'
)

image

參考來源：https://satijalab.org/signac/articles/motif_vignette.html