Motif analysis with Signac

title: "Motif analysis with Signac"

date: 'Compiled: `r format(Sys.Date(), "%B %d, %Y")`'

output: html_document

---



```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

```

In this analysis, 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](mouse_brain_vignette.html) 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:



```{r message=FALSE, warning=FALSE}

library(Signac)

library(Seurat)

library(JASPAR2020)

library(TFBSTools)

library(BSgenome.Mmusculus.UCSC.mm10)

library(patchwork)

set.seed(1234)

```



```{r message=FALSE, warning=FALSE}

mouse_brain <- readRDS("../vignette_data/adult_mouse_brain.rds")

mouse_brain

```



```{r message=FALSE, warning=FALSE}

p1 <- DimPlot(mouse_brain, label = TRUE, pt.size = 0.1) + NoLegend()

p1

```



## The Motif class



To facilitate motif analysis in Signac, we have create 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, we construct a `Motif` object and add it to our mouse brain dataset. A

motif object can be added to any Seurat assay using the `SetAssayData()`

function. See the [object interaction vignette](data_structures.html) for

more information.



```{r message=FALSE, warning=FALSE}

# Get a list of motif position frequency matrices from the JASPAR database

pfm <- getMatrixSet(

  x = JASPAR2020,

  opts = list(species = 9606, all_versions = FALSE)

)



# Scan the DNA sequence of each peak for the presence of each motif

motif.matrix <- CreateMotifMatrix(

  features = granges(mouse_brain),

  pwm = pfm,

  genome = 'mm10',

  use.counts = FALSE

)



# Create a new Mofif object to store the results

motif <- CreateMotifObject(

  data = motif.matrix,

  pwm = pfm

)



# Add the Motif object to the assay

mouse_brain <- SetAssayData(

  object = mouse_brain,

  assay = 'peaks',

  slot = 'motifs',

  new.data = motif

)

mouse_brain[["peaks"]]

```



In order to test for overrepresented motifs, we also need to compute some

sequence characteristics of the peaks, such as GC content, sequence length, and

dinucleotide frequency. The `RegionStats()` function computes this information

for us and stores the results in the feature metadata in the Seurat object.



```{r message=FALSE, warning=FALSE}

mouse_brain <- RegionStats(object = mouse_brain, genome = BSgenome.Mmusculus.UCSC.mm10)

```



## Finding overrepresented motifs



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. 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. 



```{r message=FALSE, warning=FALSE}

da_peaks <- FindMarkers(

  object = mouse_brain,

  ident.1 = 'Pvalb',

  ident.2 = 'Sst',

  only.pos = TRUE,

  test.use = 'LR',

  latent.vars = 'nCount_peaks'

)



# get top differentially accessible peaks

top.da.peak <- rownames(da_peaks[da_peaks$p_val < 0.005, ])

```



<details>

  <summary>**Optional: choosing a set of background peaks**</summary>

  

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()`.



```{r}

# 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

peaks.matched <- MatchRegionStats(

  meta.feature = GetAssayData(mouse_brain, assay = "peaks", slot = "meta.features")[open.peaks, ],

  regions = top.da.peak,

  n = 50000

)

```



`peaks.matched` can then be used as the background peak set by setting

`background=peaks.matched` in `FindMotifs()`.



</details>



```{r}

# test enrichment

enriched.motifs <- FindMotifs(

  object = mouse_brain,

  features = top.da.peak

)

```



```{r}

knitr::kable(head(enriched.motifs))

```



We can also plot the position weight matrices for the motifs, so we can

visualize the different motif sequences.



```{r}

MotifPlot(

  object = mouse_brain,

  motifs = head(rownames(enriched.motifs))

)

```



We and others have previously shown that 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 further shown that *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()`.



## Computing motif activities



We can also compute a per-cell motif activity score by running 

[chromVAR](https://greenleaflab.github.io/chromVAR/index.html). 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](https://www.nature.com/articles/nmeth.4401) for a complete description

of the method.



```{r message=FALSE, warning=FALSE, fig.width=12}

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).



```{r message=FALSE, warning=FALSE}

differential.activity <- FindMarkers(

  object = mouse_brain,

  ident.1 = 'Pvalb',

  ident.2 = 'Sst',

  only.pos = TRUE,

  test.use = 'LR',

  latent.vars = 'nCount_peaks'

)



MotifPlot(

  object = mouse_brain,

  motifs = head(rownames(differential.activity)),

  assay = 'peaks'

)

```
##################################