在本教程中，我們將學習使用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
構建Motif類
為了有助于Signac進行motif富集分析呜笑，我們需要創(chuàng)建一個Motif類來存儲所有必需的信息，其中包括位置權重矩陣（PWM）或位置頻率矩陣（PFM）的列表以及motif出現矩陣殖熟。在這里局待，我們構造了一個Motif對象，并將其添加到我們的小鼠大腦數據集中吗讶×敲停可以使用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))
)
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
我們還可以直接測試不同細胞類型之間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'
)
參考來源：https://satijalab.org/signac/articles/motif_vignette.html
