knitr::opts_chunk$set(echo = TRUE) knitr::opts_knit$set(root.dir = "../vignette_data/mouse_brain")
In this tutorial, we will perform DNA sequence motif analysis in Signac. We will explore two complementary options for performing motif analysis: one by finding overrepresented motifs in a set of differentially accessible peaks, one method performing differential motif activity analysis between groups of cells.
In this demonstration we use data from the adult mouse brain. See our vignette for the code used to generate this object, and links to the raw data. First, load the required packages and the pre-computed Seurat object:
if (!requireNamespace("TFBSTools", quietly = TRUE)) BiocManager::install("TFBSTools") if (!requireNamespace("JASPAR2020", quietly = TRUE)) BiocManager::install("JASPAR2020") if (!requireNamespace("BSgenome.Mmusculus.UCSC.mm10", quietly = TRUE)) BiocManager::install("BSgenome.Mmusculus.UCSC.mm10") if (!requireNamespace("chromVAR", quietly = TRUE)) BiocManager::install("chromVAR")
library(Signac) library(Seurat) library(JASPAR2020) library(TFBSTools) library(BSgenome.Mmusculus.UCSC.mm10) library(patchwork)
mouse_brain <- readRDS("adult_mouse_brain.rds") mouse_brain
p1 <- DimPlot(mouse_brain, label = TRUE, pt.size = 0.1) + NoLegend() p1
To add the DNA sequence motif information required for motif analyses, we can
run the AddMotifs()
function:
# Get a list of motif position frequency matrices from the JASPAR database pfm <- getMatrixSet( x = JASPAR2020, opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE) ) # add motif information mouse_brain <- AddMotifs( object = mouse_brain, genome = BSgenome.Mmusculus.UCSC.mm10, pfm = pfm )
To facilitate motif analysis in Signac, we have created the Motif
class to
store all the required information, including a list of position weight matrices
(PWMs) or position frequency matrices (PFMs) and a motif occurrence matrix.
Here, the AddMotifs()
function constructs a Motif
object and adds it to our
mouse brain dataset, along with other information such as the base composition
of each peak. A motif object can be added to any Seurat assay using the SetAssayData()
function. See the object interaction vignette for
more information.
To identify potentially important cell-type-specific regulatory sequences, we can search for DNA motifs that are overrepresented in a set of peaks that are differentially accessible between cell types.
Here, we find differentially accessible peaks between Pvalb and Sst inhibitory
interneurons. For sparse data (such as scATAC-seq), we find it is often necessary
to lower the min.pct
threshold in FindMarkers()
from the default (0.1, which
was designed for scRNA-seq data).
We then perform a hypergeometric test to test the probability of observing the motif at the given frequency by chance, comparing with a background set of peaks matched for GC content.
da_peaks <- FindMarkers( object = mouse_brain, ident.1 = 'Pvalb', ident.2 = 'Sst', only.pos = TRUE, test.use = 'LR', min.pct = 0.05, latent.vars = 'nCount_peaks' ) # get top differentially accessible peaks top.da.peak <- rownames(da_peaks[da_peaks$p_val < 0.005 & da_peaks$pct.1 > 0.2, ])
Optional: choosing a set of background peaks
Matching the set of background peaks is essential when finding enriched DNA sequence motifs. By default, we choose a set of peaks matched for GC content, but it can be sometimes be beneficial to further restrict the background peaks to those that are accessible in the groups of cells compared when finding differentially accessible peaks.
The AccessiblePeaks()
function can be used to find a set of peaks that are
open in a subset of cells. We can use this function to first restrict the set
of possible background peaks to those peaks that were open in the set of cells
compared in FindMarkers()
, and then create a GC-content-matched set of peaks
from this larger set using MatchRegionStats()
.
# find peaks open in Pvalb or Sst cells open.peaks <- AccessiblePeaks(mouse_brain, idents = c("Pvalb", "Sst")) # match the overall GC content in the peak set meta.feature <- GetAssayData(mouse_brain, assay = "peaks", layer = "meta.features") peaks.matched <- MatchRegionStats( meta.feature = meta.feature[open.peaks, ], query.feature = meta.feature[top.da.peak, ], n = 50000 )
peaks.matched
can then be used as the background peak set by setting
background=peaks.matched
in FindMotifs()
.
# test enrichment enriched.motifs <- FindMotifs( object = mouse_brain, features = top.da.peak )
knitr::kable(head(enriched.motifs))
We can also plot the position weight matrices for the motifs, so we can visualize the different motif sequences.
MotifPlot( object = mouse_brain, motifs = head(rownames(enriched.motifs)) )
Mef-family motifs, particularly
Mef2c, are enriched in Pvalb-specific peaks in scATAC-seq data
(https://doi.org/10.1016/j.cell.2019.05.031; https://doi.org/10.1101/615179),
and Mef2c is required for the development of Pvalb
interneurons (https://www.nature.com/articles/nature25999). Here our results are
consistent with these findings, and we observe a strong enrichment of Mef-family
motifs in the top results from FindMotifs()
.
We can also compute a per-cell motif activity score by running chromVAR. This allows us to visualize motif activities per cell, and also provides an alternative method of identifying differentially-active motifs between cell types.
ChromVAR identifies motifs associated with variability in chromatin accessibility between cells. See the chromVAR paper for a complete description of the method.
mouse_brain <- RunChromVAR( object = mouse_brain, genome = BSgenome.Mmusculus.UCSC.mm10 ) DefaultAssay(mouse_brain) <- 'chromvar' # look at the activity of Mef2c p2 <- FeaturePlot( object = mouse_brain, features = "MA0497.1", min.cutoff = 'q10', max.cutoff = 'q90', pt.size = 0.1 ) p1 + p2
We can also directly test for differential activity scores between cell types. This tends to give similar results as performing an enrichment test on differentially accessible peaks between the cell types (shown above).
When performing differential testing on the chromVAR z-score, we can set
mean.fxn=rowMeans
and fc.name="avg_diff"
in the FindMarkers()
function so
that the fold-change calculation computes the average difference in z-score
between the groups.
differential.activity <- FindMarkers( object = mouse_brain, ident.1 = 'Pvalb', ident.2 = 'Sst', only.pos = TRUE, mean.fxn = rowMeans, fc.name = "avg_diff" ) MotifPlot( object = mouse_brain, motifs = head(rownames(differential.activity)), assay = 'peaks' )
Session Info
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.