#' @description: Motif enrichment with chromVAR
library(Signac)
library(Seurat)
library(JASPAR2020)
library(TFBSTools)
library(BSgenome.Hsapiens.UCSC.hg38)
library(patchwork)
set.seed(101)
library(ggpubr)
library(openxlsx)
library(chromVARmotifs)
data("human_pwms_v2")
library(ComplexHeatmap)
library(circlize)
library(future)
plan("multiprocess", workers = 10)
options(future.globals.maxSize = 50000 * 1024^2) #50G
## 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.
source(file = "/home/longzhilin/Analysis_Code/Plot_colorPaletters.R")
setwd("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scATAC")
scATAC.data <- readRDS("scATAC.data.pro.rds")
Idents(scATAC.data) <- scATAC.data$AnnotatedcellType
DefaultAssay(scATAC.data) <- "Peaks"
# Get a list of motif position frequency matrices from the JASPAR database
# pfm <- getMatrixSet(
# x = JASPAR2020,
# opts = list(species = 9606, all_versions = FALSE)
# )
# add motif information
scATAC.data <- AddMotifs(
object = scATAC.data,
genome = BSgenome.Hsapiens.UCSC.hg38,
pfm = human_pwms_v2
)
# Finding overrepresented motifs
GetMotifs <- function(cellType, scATAC.object) {
print(paste0("Finding motifs for: ",cellType))
overrepresented.motifs <- FindMarkers(scATAC.object,
ident.1 = cellType,
test.use = 'LR',
min.pct = 0.05,
latent.vars = "nCount_Peaks")
enriched.motifs <- FindMotifs(object = scATAC.object, features = rownames(overrepresented.motifs[overrepresented.motifs$p_val_adj < 0.05, ]))
return(enriched.motifs)
}
# FindMarkers and write to an xlsx file with default parameters
idents <- as.character(levels(Idents(scATAC.data)))
cellType.motifs <- lapply(idents, function(x) GetMotifs(x, scATAC.object = scATAC.data))
names(cellType.motifs) <- idents
write.xlsx(cellType.motifs, file = "5.Motif/motifs.celltype.human_pwms_v2.xlsx", sheetName = idents, rowNames = T)
saveRDS(cellType.motifs, file = "5.Motif/cellType.motifs.human_pwms_v2.rds")
library(ggplot2)
library(ggseqlogo)
PWMatrixToProbMatrix <- function(x){
if (class(x) != "PWMatrix") stop("x must be a TFBSTools::PWMatrix object")
m <- (exp(as(x, "matrix"))) * TFBSTools::bg(x)/sum(TFBSTools::bg(x))
m <- t(t(m)/colSums(m))
m
}
pdf("5.Motif/cellType.specific.TFs.logo.pdf")
cellType.TFs.PPM <- sapply(idents, function(x){
cellType.top.TFs <- head(rownames(cellType.motifs[[x]]))
cellType.top.TFs <- gsub("-", "_", cellType.top.TFs)
index <- match(cellType.top.TFs, names(human_pwms_v2))
PPM.list <- lapply(index, function(y){
PPM <- PWMatrixToProbMatrix(human_pwms_v2[[y]])
})
names(PPM.list) <- as.character(unlist(scATAC.data@assays$Peaks@motifs@motif.names[index]))
p <- ggseqlogo(PPM.list) + ggtitle(x) + theme(plot.title = element_text(hjust = 0.5))
print(p)
return(PPM.list)
})
dev.off()
##########################################ChromVAR
#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.
# Computing motif activities
scATAC.data <- RunChromVAR(
object = scATAC.data,
genome = BSgenome.Hsapiens.UCSC.hg38
)
DefaultAssay(scATAC.data) <- 'chromvar'
#cellType.chromvar.activity <- FindAllMarkers(scATAC.data, group.by = "AnnotatedcellType", test.use = 'LR', latent.vars = "nCount_Peaks")
library(tidyverse)
GetChromvarActivities <- function(cellType, scATAC.object, motif.info) {
print(paste0("Finding chromVAR activities for: ",cellType))
differential.activity <- FindMarkers(scATAC.object,
ident.1 = cellType,
test.use = 'LR',
logfc.threshold = 0,
latent.vars = "nCount_Peaks")
motifs <- gsub("-", "_", rownames(differential.activity))
motifNames <- sapply(motifs, function(x) motif.info@motif.names[[x]])
return(cbind(differential.activity, gene = motifNames))
}
idents <- as.character(levels(Idents(scATAC.data)))
cellType.motifs.chromVAR <- lapply(idents, function(x) GetChromvarActivities(x, scATAC.object = scATAC.data, motif.info = scATAC.data@assays$Peaks@motifs))
names(cellType.motifs.chromVAR) <- idents
cellType.motifs.chromVAR <- lapply(cellType.motifs.chromVAR, function(x){
up.x <- x %>% filter(avg_log2FC>0) %>% arrange(desc(avg_log2FC))
down.x <- x %>% filter(avg_log2FC<=0) %>% arrange(avg_log2FC)
x <- rbind(up.x, down.x)
return(x)
})
names(cellType.motifs.chromVAR) <- idents
write.xlsx(cellType.motifs.chromVAR, file = "5.Motif/cellType.motifs.chromVAR.human_pwms_v2.xlsx", sheetName = idents, rowNames = T)
saveRDS(cellType.motifs.chromVAR, file = "5.Motif/cellType.motifs.chromVAR.human_pwms_v2.rds")
pdf("5.Motif/cellType.specific.TFs.chromVAR.logo.pdf")
cellType.TFs.PPM <- sapply(idents, function(x){
cellType.top.TFs <- head(rownames(cellType.motifs.chromVAR[[x]]))
cellType.top.TFs <- gsub("-", "_", cellType.top.TFs)
index <- match(cellType.top.TFs, names(human_pwms_v2))
PPM.list <- lapply(index, function(y){
PPM <- PWMatrixToProbMatrix(human_pwms_v2[[y]])
})
names(PPM.list) <- as.character(unlist(scATAC.data@assays$Peaks@motifs@motif.names[index]))
p <- ggseqlogo(PPM.list) + ggtitle(x) + theme(plot.title = element_text(hjust = 0.5))
print(p)
return(PPM.list)
})
dev.off()
saveRDS(scATAC.data, "scATAC.data.pro.rds")
Datasets
Models
