Skip to main content

SCmultiomics Analysis in LongCOVID

21 Motif Enrichment

Once peaks are linked to genes, transcription factor (TF) binding motifs are scanned across all ATAC peaks using the JASPAR2020 core vertebrate motif database. Motif presence per peak enables two downstream analyses: (1) per-cluster motif enrichment testing (this chapter), which asks which TF motifs are over-represented in a cluster’s accessible peaks relative to a matched background, and (2) chromVAR TF activity scoring (next chapter), which estimates per-cell TF activity directly from motif accessibility.

21.1 Load peak-linked object 

library(Seurat)
library(Signac)
library(TFBSTools)
library(JASPAR2020)
library(BSgenome.Hsapiens.UCSC.hg38)
library(dplyr)


mygenome <- BSgenome.Hsapiens.UCSC.hg38
genome(mygenome) <- NA

21.2 Add motif matches

AddMotifs() scans every peak in the ATAC assay for matches to each motif in the JASPAR CORE vertebrate collection and stores a peak-by-motif match matrix in the assay. A lookup table mapping JASPAR motif IDs to gene symbols is built alongside it, since AddMotifs() output is indexed by motif ID, not TF symbol.

DefaultAssay(combined) <- "ATAC"

pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(collection = "CORE", tax_group = "vertebrates", all_versions = FALSE)
)

df_pfm <- data.frame(t(sapply(pfm, function(x)
  c(id = x@ID, name = x@name, symbol = ifelse(!is.null(x@tags$symbol), x@tags$symbol, NA)))))

combined <- AddMotifs(combined, genome = mygenome, pfm = pfm)

21.3 Motif enrichment per cluster

Enrichment is tested against a background of peaks matched on GC content, peak width, and accessibility rank (MatchRegionStats()), so that enrichment reflects motif over-representation rather than generic peak-composition differences. Motif enrichment is tested across all peaks used for cluster markers.

DefaultAssay(combined) <- "ATAC"

open_peaks <- AccessiblePeaks(combined)

peaks_matched <- MatchRegionStats(meta.feature = combined[['ATAC']]@meta.features[open_peaks, ],
                                  query.feature = combined[['ATAC']]@meta.features,
                                  n = 50000)

motif_enrichment <- FindMotifs(combined,
                               features = rownames(combined),
                               background = peaks_matched) %>%
  mutate(symbol = setNames(ifelse(is.na(df_pfm$symbol), df_pfm$name, df_pfm$symbol), df_pfm$id)[motif])

21.4 Top enriched motifs 

MotifPlot(combined, motifs = motif_enrichment$motif[1:10], ncol = 4)

Download PDF

21.5 Save motif-annotated object 

SaveSeuratRds(combined, "data/seurat_object_cross_modal_integration_peakslinks_motifs.rds")

# FIXED: original write.table() call used a relative path ("motifs.txt"); output
# paths should be absolute and under the shared processing directory.
write.table(motif_enrichment,
           "data/motifs.txt",
           row.names = FALSE)