Question about Computing motif activities

[Sophiesze]

Hi Tim,

I'm trying to identify the motif difference between patients and healthy controls in CD4+ cells by using RunChromvar. My code is attached below. It is strange that although the avg_logFC seems high, the Heatmap and VlnPlot of the motif doesn't look like have any differences. I hope you can help me with the problem I met. Thanks!

Best, Sophie

####
pwm <- getMatrixSet(
  x = JASPAR2018,
  opts = list(species = 9606, all_versions = FALSE)
)
#### motif
seurat=subset(seurat,features=VariableFeatures(seurat))
motif.matrix <- CreateMotifMatrix(
  features = StringToGRanges(rownames(seurat), sep = c(":", "-")),
  pwm = pwm,
  genome = BSgenome.Hsapiens.UCSC.hg38,
  sep = c(":", "-")
)

# Create a new Mofif object to store the results
motif <- CreateMotifObject(
  data = motif.matrix,
  pwm = pwm
)
# Add the Motif object to the assay
seurat[['peaks']] <- AddMotifObject(
  object = seurat[['peaks']],
  motif.object = motif
)
seurat <- RegionStats(
  object = seurat,
  genome = BSgenome.Hsapiens.UCSC.hg38,
  sep = c(":", "-")
)

#Computing with ChromVAR
print('Computing with ChromVAR')
seurat <- RunChromVAR(
  object = seurat,
  genome = BSgenome.Hsapiens.UCSC.hg38
)

DefaultAssay(seurat) <- 'chromvar'
differential.activity <- FindAllMarkers(seurat,only.pos=FALSE, min.pct= 0.25,logfc.threshold=0.25,test.use='LR',latent.vars='nCount_peaks')

seurat <- ScaleData(seurat)
top10 <- differential.activity %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(seurat, features = top10$gene) + NoLegend()

seurat
An object of class Seurat 
101140 features across 6614 samples within 4 assays 
Active assay: new_chromvar (452 features, 0 variable features)
 3 other assays present: peaks, RNA, chromvar
 4 dimensional reductions calculated: lsi, umap, tsne, harmony
####The assay 'new_chromvar' is the same as the assay 'chromvar' with different rownames(with motif.names instead of Matrix ID).

[timoast]

How did you set the cell identities in the object before doing the differential expression test?

[Sophiesze]

I defined CD4+ cells only by promoters and motifs from previous studies done by CHIP-seq or bulk ATAC seq, and so do other cell types. I did not integrate scATAC data with sc-RNA data.

[timoast]

I meant more how are you setting the identities in the object? My guess is that you're not comparing the right groups of cells in the violin plot / heatmap and in the FindMarkers function. In the heatmap you can see within the red/blue groups (name is cropped out) there are different populations of cells that are high/low for FOS::JUND. Can you post the full you that you ran?

[Sophiesze]

Codes are attached. Thanks!

CODE.zip

[timoast]

It's a bit hard to follow all the code because it's in different files and there are a lot of things going on, in general it's helpful to provide a minimal reproducible example.

My understanding is that you identified CD4+ T-cells, subset the object to only contain those cells, and set the identities in the subset object to either "HC" or "PATIENT" depending on some additional metadata. You then found differentially active motifs by running FindMarkers on the chromVAR assay.

I would try to verify again that you have the identities set correctly when running VlnPlot and DoHeatmap, as it looks like there are separate populations of cells that are high and low for this motif, which is more consistent with the fold change reported in the DE test. Also, for the heatmap you probably want to use the data rather than scale.data slot when plotting the chromVAR assay.

[Sophiesze]

Thanks for your advice!