0.1 Set up

We have some modified function which is provided in our github repository(https://github.com/liranmao/Spatial_multi_omics/tree/main) and you have to load them first.

rm(list=ls())
library(ArchR)
library(Seurat)
library(grid)
library(BSgenome.Mmusculus.UCSC.mm10)
library(ggplot2)
library(patchwork)
library(dplyr)
library(gridExtra)
library(pander)
# options("Seurat.object.assay.version" = "v3")
# These files can be found in Github repo
source('./scripts/getGeneScore_ArchR.R')
source('./scripts/SpatialPlot_new.R')

0.2 Process the fragment file for each modality

First we should process each modality seperately. Here we stated we H3K27me3 using ArchR package. Then create folder for each modality and put the spatial and fragment file inside each folder. For this demo, we have provided all the input files here for easier exploration:.

########## archr project creation
threads = 8
addArchRThreads(threads = threads)
addArchRGenome("mm10")

setwd('/mnt/HDD1/Users/liran/02_nano_review/17nanobody_deep/generate_RMD/Spananob17_H3K27ac')
for (i in c('markers_list', 'all_plot')){
  folder_path <- paste0("./",i)
  dir.create(folder_path)
}

sampleNames <- 'Spananob17_H3K27ac'
inputFiles <- './Spnanob17.fragments.sort.bed.gz'


ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = sampleNames,
  filterTSS = 0, #Dont set this too high because you can always increase later
  filterFrags = 0,
  minFrags = 0,
  maxFrags = 1e+07,
  addTileMat = TRUE,
  addGeneScoreMat = TRUE,
  offsetPlus = 0,
  offsetMinus = 0,
  TileMatParams = list(tileSize = 5000)
)
ArrowFiles

projHeme1 <- ArchRProject(
  ArrowFiles = ArrowFiles, 
  outputDirectory = sampleNames,
  copyArrows = TRUE #This is recommened so that if you modify the Arrow files you have an original copy for later usage.
)
projHeme1

saveArchRProject(ArchRProj = projHeme1, outputDirectory = paste0("Save-", sampleNames), load = FALSE)
projHeme1 <- loadArchRProject(path = paste0("Save-", sampleNames), force = FALSE, showLogo = TRUE)


############ Combine with spatial information
## Prepare meta data
meta.data <- as.data.frame(getCellColData(ArchRProj = projHeme1))
meta.data['cellID_archr'] <- row.names(meta.data)
new_row_names <- row.names(meta.data)
new_row_names <- unlist(lapply(new_row_names, function(x) gsub(".*#","", x)))
new_row_names <- unlist(lapply(new_row_names, function(x) gsub("-.*","", x)))
row.names(meta.data) <- new_row_names

data.dir <- getwd()
assay = "Spatial"
filter.matrix = TRUE
slice = "slice1"

# There is a spatial image for each slide
image <- Read10X_Image(image.dir = file.path(data.dir, "spatial"), 
                       filter.matrix = filter.matrix)

## export meta data
meta.data.spatial <- meta.data[row.names(image@coordinates), ]

write.table(meta.data.spatial, paste0(sampleNames, '_meta_data.txt'), row.names = TRUE, col.names = TRUE)

## filter off-tissue tixels using image data
projCUTA <- projHeme1[meta.data.spatial$cellID_archr, ]
projCUTA

p <- plotFragmentSizes(ArchRProj = projCUTA)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_FragmentSizes.pdf'), width = 5, height = 5)
print(p)
dev.off()

########### dimension reduction, clustering, and add UMAP embedding
projCUTA <- addIterativeLSI(
  ArchRProj = projCUTA,
  useMatrix = "TileMatrix", 
  name = "IterativeLSI", 
  iterations = 2, 
  clusterParams = list( 
    resolution = c(0.2), 
    sampleCells = 10000, 
    n.start = 10
  ), 
  varFeatures = 25000, 
  dimsToUse = 1:30,
  force = TRUE
)

projCUTA <- addClusters(
  input = projCUTA,
  reducedDims = "IterativeLSI",
  method = "Seurat",
  name = "Clusters",
  resolution = 1,
  force = TRUE
)

projCUTA <- addUMAP(
  ArchRProj = projCUTA, 
  reducedDims = "IterativeLSI", 
  name = "UMAP", 
  nNeighbors = 30, 
  minDist = 0.5, 
  metric = "cosine",
  force = TRUE
)

projCUTA <- addImputeWeights(projCUTA)

saveArchRProject(ArchRProj = projCUTA, outputDirectory = paste0("Save-inTissue-", sampleNames), load = FALSE)
projCUTA <- loadArchRProject(path = paste0("Save-inTissue-", sampleNames), force = FALSE, showLogo = TRUE)
projCUTA

## QC plot
df <- getCellColData(projCUTA, select = c("log10(nFrags)", "TSSEnrichment"))
p1 <- ggPoint(
  x = df[,1], 
  y = df[,2], 
  colorDensity = TRUE,
  continuousSet = "sambaNight",
  xlabel = "Log10 Unique Fragments",
  ylabel = "TSS Enrichment",
  xlim = c(log10(100), 4.5),
  ylim = c(1, 8),
  baseSize = 12
)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, "TSS-vs-Frags.pdf"))
print(p1)
dev.off()

markersGS <- getMarkerFeatures(
  ArchRProj = projCUTA, 
  useMatrix = "GeneScoreMatrix", 
  groupBy = "Clusters",
  testMethod = "wilcoxon"
)

markerList_pos <- getMarkers(markersGS, cutOff = "FDR < 0.05 & Log2FC >= 0.5")

for (i in seq_len(length(markerList_pos))) {
  write.table(markerList_pos[[i]], file=paste0('./markers_list/', sampleNames, '_C', i, '_markers.txt'),
              quote=FALSE, sep='\t', col.names = TRUE, row.names = TRUE)
}

heatmapGS <- plotMarkerHeatmap(
  seMarker = markersGS,
  cutOff = "FDR <= 0.05 & Log2FC <= -1",
  transpose = TRUE
)

heatmapGS_complex <- ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")

# Save the heatmap as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, "matker_motif_heatmap.pdf"))
print(heatmapGS_complex)
dev.off()

markerGenes <- list()
for (i in seq_len(length(markerList_pos))) {
  markerGenes <- c(markerGenes, markerList_pos[[i]]$name)
}
markerGenes <- unlist(markerGenes)


############ Spatial plots
meta.data <- as.data.frame(getCellColData(ArchRProj = projCUTA))
meta.data['cellID_archr'] <- row.names(meta.data)
new_row_names <- row.names(meta.data)
new_row_names <- unlist(lapply(new_row_names, function(x) gsub(".*#","", x)))
new_row_names <- unlist(lapply(new_row_names, function(x) gsub("-.*","", x)))
row.names(meta.data) <- new_row_names

projCUTA <- addImputeWeights(projCUTA)
gene_score <- getGeneScore_ArchR(ArchRProj = projCUTA, name = markerGenes, imputeWeights = getImputeWeights(projCUTA))
saveRDS(gene_score, paste0(sampleNames,'_gene_score.rds'))
gene_score <- readRDS(paste0(sampleNames,'_gene_score.rds'))

## create seurat object
object <- CreateSeuratObject(counts = gene_score, assay = assay, meta.data = meta.data)

image <- Read10X_Image(image.dir = file.path(data.dir, "spatial"), filter.matrix = filter.matrix)
image <- image[Cells(x = object)]
DefaultAssay(object = image) <- assay
object[[slice]] <- image

spatial.obj <- object

p2 <- VlnPlot(spatial.obj, features = "nFrags", pt.size = 0.1, log = TRUE) + NoLegend()
median(meta.data$nFrags)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_nfrags.pdf'), width = 5, height = 5)
print(p2)
dev.off()

p3 <- SpatialPlot_new(spatial.obj, features = "nFrags", pt.size.factor = 4.5, min.cutoff = "q10", max.cutoff = "q90", image.alpha = 0, stroke = 0) + 
  theme(legend.position = "right", legend.text=element_text(size=15), legend.title=element_text(size=15))
p3$layers[[1]]$aes_params <- c(p3$layers[[1]]$aes_params, shape=22)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_nFrags_spatial.pdf'), width = 5, height = 5)
print(p3)
dev.off()

n_clusters <- length(unique(projCUTA$Clusters))
cols <- ArchRPalettes$stallion[as.character(seq_len(n_clusters))]
names(cols) <- paste0('C', seq_len(n_clusters))
cols

p4 <- SpatialPlot(spatial.obj, label = FALSE, label.size = 3, group.by = 'Clusters', pt.size.factor = 4.5, cols = cols, image.alpha = 0, stroke = 0)
p4$layers[[1]]$aes_params <- c(p4$layers[[1]]$aes_params, shape=22)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_clusters_spatial.pdf'), width = 5, height = 5)
print(p4)
dev.off()

p5 <- plotEmbedding(ArchRProj = projCUTA, colorBy = "cellColData", name = "Clusters", embedding = "UMAP", size = 0.5)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_clusters_umap.pdf'), width = 5, height = 5)
print(p5)
dev.off()

Do the same for H3K27me3. You only have to change the working directory, sample names and input files. Run the first chunk again to clear the previous workspace.

########## archr project creation
threads = 8
addArchRThreads(threads = threads)
addArchRGenome("mm10")

setwd('/mnt/HDD1/Users/liran/02_nano_review/17nanobody_deep/generate_RMD/Spananob17_H3K27me3')
for (i in c('markers_list', 'all_plot')){
  folder_path <- paste0("./",i)
  dir.create(folder_path)
}

sampleNames <- 'Spananob17_H3K27me3'
inputFiles <- './Spnanob17.fragments.sort.bed.gz'


ArrowFiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = sampleNames,
  filterTSS = 0, #Dont set this too high because you can always increase later
  filterFrags = 0,
  minFrags = 0,
  maxFrags = 1e+07,
  addTileMat = TRUE,
  addGeneScoreMat = TRUE,
  offsetPlus = 0,
  offsetMinus = 0,
  TileMatParams = list(tileSize = 5000)
)
ArrowFiles

projHeme1 <- ArchRProject(
  ArrowFiles = ArrowFiles, 
  outputDirectory = sampleNames,
  copyArrows = TRUE #This is recommened so that if you modify the Arrow files you have an original copy for later usage.
)
projHeme1

saveArchRProject(ArchRProj = projHeme1, outputDirectory = paste0("Save-", sampleNames), load = FALSE)
projHeme1 <- loadArchRProject(path = paste0("Save-", sampleNames), force = FALSE, showLogo = TRUE)


############ Combine with spatial information
## Prepare meta data
meta.data <- as.data.frame(getCellColData(ArchRProj = projHeme1))
meta.data['cellID_archr'] <- row.names(meta.data)
new_row_names <- row.names(meta.data)
new_row_names <- unlist(lapply(new_row_names, function(x) gsub(".*#","", x)))
new_row_names <- unlist(lapply(new_row_names, function(x) gsub("-.*","", x)))
row.names(meta.data) <- new_row_names

data.dir <- getwd()
assay = "Spatial"
filter.matrix = TRUE
slice = "slice1"

# There is a spatial image for each slide
image <- Read10X_Image(image.dir = file.path(data.dir, "spatial"), 
                       filter.matrix = filter.matrix)

## export meta data
meta.data.spatial <- meta.data[row.names(image@coordinates), ]

write.table(meta.data.spatial, paste0(sampleNames, '_meta_data.txt'), row.names = TRUE, col.names = TRUE)

## filter off-tissue tixels using image data
projCUTA <- projHeme1[meta.data.spatial$cellID_archr, ]
projCUTA

p <- plotFragmentSizes(ArchRProj = projCUTA)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_FragmentSizes.pdf'), width = 5, height = 5)
print(p)
dev.off()

########### dimension reduction, clustering, and add UMAP embedding
projCUTA <- addIterativeLSI(
  ArchRProj = projCUTA,
  useMatrix = "TileMatrix", 
  name = "IterativeLSI", 
  iterations = 2, 
  clusterParams = list( 
    resolution = c(0.2), 
    sampleCells = 10000, 
    n.start = 10
  ), 
  varFeatures = 25000, 
  dimsToUse = 1:30,
  force = TRUE
)

projCUTA <- addClusters(
  input = projCUTA,
  reducedDims = "IterativeLSI",
  method = "Seurat",
  name = "Clusters",
  resolution = 1,
  force = TRUE
)

projCUTA <- addUMAP(
  ArchRProj = projCUTA, 
  reducedDims = "IterativeLSI", 
  name = "UMAP", 
  nNeighbors = 30, 
  minDist = 0.5, 
  metric = "cosine",
  force = TRUE
)

projCUTA <- addImputeWeights(projCUTA)

saveArchRProject(ArchRProj = projCUTA, outputDirectory = paste0("Save-inTissue-", sampleNames), load = FALSE)
projCUTA <- loadArchRProject(path = paste0("Save-inTissue-", sampleNames), force = FALSE, showLogo = TRUE)
projCUTA

## QC plot
df <- getCellColData(projCUTA, select = c("log10(nFrags)", "TSSEnrichment"))
p1 <- ggPoint(
  x = df[,1], 
  y = df[,2], 
  colorDensity = TRUE,
  continuousSet = "sambaNight",
  xlabel = "Log10 Unique Fragments",
  ylabel = "TSS Enrichment",
  xlim = c(log10(100), 4.5),
  ylim = c(1, 8),
  baseSize = 12
)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, "TSS-vs-Frags.pdf"))
print(p1)
dev.off()

markersGS <- getMarkerFeatures(
  ArchRProj = projCUTA, 
  useMatrix = "GeneScoreMatrix", 
  groupBy = "Clusters",
  testMethod = "wilcoxon"
)

markerList_pos <- getMarkers(markersGS, cutOff = "FDR < 0.05 & Log2FC >= 0.5")

for (i in seq_len(length(markerList_pos))) {
  write.table(markerList_pos[[i]], file=paste0('./markers_list/', sampleNames, '_C', i, '_markers.txt'),
              quote=FALSE, sep='\t', col.names = TRUE, row.names = TRUE)
}

heatmapGS <- plotMarkerHeatmap(
  seMarker = markersGS,
  cutOff = "FDR <= 0.05 & Log2FC <= -1",
  transpose = TRUE
)

heatmapGS_complex <- ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")

# Save the heatmap as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, "matker_motif_heatmap.pdf"))
print(heatmapGS_complex)
dev.off()

markerGenes <- list()
for (i in seq_len(length(markerList_pos))) {
  markerGenes <- c(markerGenes, markerList_pos[[i]]$name)
}
markerGenes <- unlist(markerGenes)


############ Spatial plots
meta.data <- as.data.frame(getCellColData(ArchRProj = projCUTA))
meta.data['cellID_archr'] <- row.names(meta.data)
new_row_names <- row.names(meta.data)
new_row_names <- unlist(lapply(new_row_names, function(x) gsub(".*#","", x)))
new_row_names <- unlist(lapply(new_row_names, function(x) gsub("-.*","", x)))
row.names(meta.data) <- new_row_names

projCUTA <- addImputeWeights(projCUTA)
gene_score <- getGeneScore_ArchR(ArchRProj = projCUTA, name = markerGenes, imputeWeights = getImputeWeights(projCUTA))
saveRDS(gene_score, paste0(sampleNames,'_gene_score.rds'))
gene_score <- readRDS(paste0(sampleNames,'_gene_score.rds'))

## create seurat object
object <- CreateSeuratObject(counts = gene_score, assay = assay, meta.data = meta.data)

image <- Read10X_Image(image.dir = file.path(data.dir, "spatial"), filter.matrix = filter.matrix)
image <- image[Cells(x = object)]
DefaultAssay(object = image) <- assay
object[[slice]] <- image

spatial.obj <- object

p2 <- VlnPlot(spatial.obj, features = "nFrags", pt.size = 0.1, log = TRUE) + NoLegend()
median(meta.data$nFrags)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_nfrags.pdf'), width = 5, height = 5)
print(p2)
dev.off()

p3 <- SpatialPlot_new(spatial.obj, features = "nFrags", pt.size.factor = 4.5, min.cutoff = "q10", max.cutoff = "q90", image.alpha = 0, stroke = 0) + 
  theme(legend.position = "right", legend.text=element_text(size=15), legend.title=element_text(size=15))
p3$layers[[1]]$aes_params <- c(p3$layers[[1]]$aes_params, shape=22)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_nFrags_spatial.pdf'), width = 5, height = 5)
print(p3)
dev.off()

n_clusters <- length(unique(projCUTA$Clusters))
cols <- ArchRPalettes$stallion[as.character(seq_len(n_clusters))]
names(cols) <- paste0('C', seq_len(n_clusters))
cols

p4 <- SpatialPlot(spatial.obj, label = FALSE, label.size = 3, group.by = 'Clusters', pt.size.factor = 4.5, cols = cols, image.alpha = 0, stroke = 0)
p4$layers[[1]]$aes_params <- c(p4$layers[[1]]$aes_params, shape=22)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_clusters_spatial.pdf'), width = 5, height = 5)
print(p4)
dev.off()

p5 <- plotEmbedding(ArchRProj = projCUTA, colorBy = "cellColData", name = "Clusters", embedding = "UMAP", size = 0.5)

# Save the plot as PDF in 'all_plot'
pdf(file = paste0("all_plot/", sampleNames, '_clusters_umap.pdf'), width = 5, height = 5)
print(p5)
dev.off()

To reproduce exact the same results from our paper, we have provided our archR object and you can directly load them through the following code.

path <- pathe_to_archR_object, eg. '/mnt/HDD1/Users/liran/02_nano_review/17nanobody_gp/Save-inTissue-Spananob17_H3K27ac/'
projCUTA <- loadArchRProject(path = path, force = FALSE, showLogo = TRUE)
projCUTA

0.3 Multi-omics integration

For our multi-omics data integration, we consolidated all modalities into a single Seurat object. The ATAC and CUT&Tag data integration utilized a 500bp peak matrix generated by addReproduciblePeakSet from ArchR, applying Macs2 for peak calling. RNA data integration was based on a log-normalized gene expression matrix. We applied Weighted Nearest Neighbors (WNN) analysis with FindMultiModalNeighbors for clustering, utilizing UMAP and spatial mapping for visualization. Subsequently, cell type clusters were refined through FindSubCluster within Seurat, based on the wsnn graph. This streamlined approach facilitated a precise analysis of cellular heterogeneity within the multi-omics dataset.

############## Basic settings
threads = 8
addArchRThreads(threads = threads)
addArchRGenome("mm10")

samplepath_list <- c('/mnt/HDD1/Users/liran/02_nano_review/17nanobody_deep_ori/Spnanob17_h3k27me3_deep/Save-inTissue-Spananob17_H3K27me3/',
                     '/mnt/HDD1/Users/liran/02_nano_review/17nanobody_deep_ori/Spnanob17_h3k27ac_deep/Save-inTissue-Spananob17_H3K27ac/')
modalitylist <- c('BC6_h3k27me3', 'BC7_h3k27ac')
colnameslist<-c('Spananob17_H3K27me3#','Spananob17_H3K27ac#')

for (i in c('combine_plot')){
  folder_path <- paste0("./",i)
  # Create the new folder
  dir.create(folder_path)
}

################ Call peaks through archR and create the seurat object
for (i in 1:2){
  samplepath <- samplepath_list[i]
  modality <- modalitylist[i]
  colname <- colnameslist[i]
  projCUTA <- loadArchRProject(path = samplepath, force = FALSE, showLogo = TRUE)

  projCUTA <- addGroupCoverages(ArchRProj = projCUTA, groupBy = "Clusters")
  pathToMacs2 <- '/home/liran/miniconda3/envs/macs2/bin/macs2'
  projCUTA <- addReproduciblePeakSet(
    ArchRProj = projCUTA, 
    groupBy = "Clusters", 
    pathToMacs2 = pathToMacs2,
    force = TRUE
  )
  
  projCUTA <- addPeakMatrix(projCUTA)

  PeakSet <- getPeakSet(projCUTA)
  SummarizedExperiment<- getMatrixFromProject( ArchRProj = projCUTA, useMatrix = "PeakMatrix")
  SummarizedExperiment_assay <- assay(SummarizedExperiment)
  
  # Combine seqnames and ranges into column names
  col_names <- paste(seqnames(PeakSet), ranges(PeakSet), sep = "-")
  # Set the column names of the sparse matrix
  rownames(SummarizedExperiment_assay) <- col_names
  
  colnames_sparse <- colnames(SummarizedExperiment_assay)
  new_colnames <- sub(colname, "", colnames_sparse)
  new_colnames <- sub("-1", "", new_colnames)
  colnames(SummarizedExperiment_assay) <-new_colnames
  if(i == 1){
    seurat.wnn.peak <- CreateSeuratObject(counts = SummarizedExperiment_assay, assay = paste0('peaks_',modality))
  } else {
    # create a new assay to store ADT information
    adt_assay <- CreateAssayObject(counts = SummarizedExperiment_assay)
    adt_assay <- subset(adt_assay, cells = colnames(seurat.wnn.peak[[paste0('peaks_',modalitylist[1])]]))
    # add this assay to the previously created Seurat object
    seurat.wnn.peak[[paste0('peaks_',modality)]] <- adt_assay
  }
  
}

########### Multi-omics integration through WNN
for(x in modalitylist){
  DefaultAssay(seurat.wnn.peak) <- paste0('peaks_',x)
  seurat.wnn.peak <- RunTFIDF(seurat.wnn.peak) %>% FindTopFeatures() %>% RunSVD(reduction.name = paste0(x,'_lsi'))
}

seurat.wnn.peak <- FindMultiModalNeighbors(
  seurat.wnn.peak, reduction.list = list( 'BC6_h3k27me3_lsi', 'BC7_h3k27ac_lsi'), 
  dims.list = list(2:30,2:30), modality.weight.name = "histone.weight",k.nn = 10
)

seurat.wnn.peak <- FindClusters(seurat.wnn.peak, graph.name = "wsnn", algorithm = 3, resolution = 1, verbose = FALSE)


for(x in modalities){
  assay <- paste0('peaks_',x,'.weight')
  p2 <- FeaturePlot(seurat.wnn.peak,assay,min.cutoff = 0.2,max.cutoff = 0.6) + scale_color_viridis_c()
  pdf(paste0("combine_plot/", outputname, '_', x,"_clusters_umap_wsnn_weight.pdf"), width = 5, height = 5)
  print(p2)
  dev.off()
}

n_clusters <- length(unique(seurat.wnn.peak$wsnn_res.1))
cols <- ArchRPalettes$stallion[as.character(seq_len(n_clusters))]
names(cols) <- paste0( seq_len(n_clusters)-1)
cols

p3 <- DimPlot(seurat.wnn.peak,label=TRUE,pt.size=0.5, cols = cols)

pdf(paste0("combine_plot/", outputname, "_clusters_umap_wsnn_res.1.pdf"), width = 5.3, height = 4)
print(p3)
dev.off()

p3

You can also visualize the spatial plot.

## add image, use any path to spatial folder
data.dir <- "/mnt/HDD1/Users/liran/02_nano_review/17nanobody_deep/generate_RMD/Spananob17_H3K27me3/"
filter.matrix = TRUE
image <- Read10X_Image(image.dir = file.path(data.dir, "spatial"), filter.matrix = filter.matrix)
image <- image[sub("-1", "", sub("Spananob18_BC5#", "", Cells(x = seurat.wnn.peak)))]
DefaultAssay(seurat.wnn.peak = image) <- 'peaks_BC7_h3k27ac'
slice = "slice1"
seurat.wnn.peak[[slice]] <- image


# spatial plot weight cluster
for(x in modalities){
  assay <- paste0('peaks_',x,'.weight')
  
  p3 <- SpatialPlot(seurat.wnn.peak, label = FALSE, label.size = 3, features = assay,pt.size.factor = 4.5,  image.alpha = 0, stroke = 0, min.cutoff = 0.3, max.cutoff = 0.9)+ scale_color_viridis_c()
  p3$layers[[1]]$aes_params <- c(p3$layers[[1]]$aes_params, shape=22)
  
  pdf(paste0("combine_plot/", outputname, '_', x,"_clusters_spatial_wsnn_weight.pdf"), width = 5, height = 5)
  print(p3)
  dev.off()
}


# spatial plot cluster
p4 <- SpatialPlot(seurat.wnn.peak, label = FALSE, label.size = 3, group.by = 'wsnn_res.1', pt.size.factor = 4.5, cols = cols, image.alpha = 0, stroke = 0)
p4$layers[[1]]$aes_params <- c(p4$layers[[1]]$aes_params, shape=22)
p4

pdf(paste0("combine_plot/", outputname, "_clusters_spatial_wsnn_res.1.pdf"), width = 5, height = 5)
print(p4)
dev.off()

# Save the Seurat object to an RDS file
saveRDS(seurat.wnn.peak, file = "./seurat.wnn.peak.rds")

p4

To reproduce exact the same results from our paper, we have provided our Seurat object with all omics layers and you can directly load them through the following code.

spatial.obj <- readRDS('/mnt/HDD1/Users/liran/02_nano_review/17nanobody_deep/combine_17nanobody_deep/seurat.wnn.peak_upload.rds')

p1 <- DimPlot(spatial.obj, reduction = "wnn.umap", label = TRUE)
p1

n_clusters <- length(unique(spatial.obj$wsnn_res.1))
cols <- ArchRPalettes$stallion[as.character(seq_len(n_clusters))]
names(cols) <- paste0( seq_len(n_clusters)-1)

p1 <- SpatialPlot(spatial.obj, label = FALSE, label.size = 3, group.by = 'wsnn_res.1', pt.size.factor = 5, cols = cols,  image.alpha = 0, stroke = 0)
p1$layers[[1]]$aes_params <- c(p1$layers[[1]]$aes_params, shape=22)
p1

0.4 Session Info

pander(sessionInfo())

R version 4.4.1 (2024-06-14)

Platform: x86_64-pc-linux-gnu

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats4, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: pander(v.0.6.5), dplyr(v.1.1.4), patchwork(v.1.2.0), BSgenome.Mmusculus.UCSC.mm10(v.1.4.3), BSgenome(v.1.72.0), rtracklayer(v.1.64.0), BiocIO(v.1.14.0), Biostrings(v.2.72.1), XVector(v.0.44.0), SeuratObject(v.4.1.4), Seurat(v.4.4.0), rhdf5(v.2.48.0), SummarizedExperiment(v.1.34.0), Biobase(v.2.64.0), MatrixGenerics(v.1.16.0), Rcpp(v.1.0.12), Matrix(v.1.7-0), GenomicRanges(v.1.56.1), GenomeInfoDb(v.1.40.1), IRanges(v.2.38.0), S4Vectors(v.0.42.0), BiocGenerics(v.0.50.0), matrixStats(v.1.3.0), data.table(v.1.15.4), stringr(v.1.5.1), plyr(v.1.8.9), magrittr(v.2.0.3), ggplot2(v.3.5.1), gtable(v.0.3.5), gtools(v.3.9.5), gridExtra(v.2.3) and ArchR(v.1.0.2)

loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), rstudioapi(v.0.16.0), jsonlite(v.1.8.8), spatstat.utils(v.3.0-5), rmarkdown(v.2.27), zlibbioc(v.1.50.0), vctrs(v.0.6.5), ROCR(v.1.0-11), Rsamtools(v.2.20.0), spatstat.explore(v.3.2-7), RCurl(v.1.98-1.14), htmltools(v.0.5.8.1), S4Arrays(v.1.4.1), curl(v.5.2.1), Rhdf5lib(v.1.26.0), SparseArray(v.1.4.8), sass(v.0.4.9), sctransform(v.0.4.1), parallelly(v.1.37.1), KernSmooth(v.2.23-24), bslib(v.0.7.0), htmlwidgets(v.1.6.4), ica(v.1.0-3), plotly(v.4.10.4), zoo(v.1.8-12), cachem(v.1.1.0), GenomicAlignments(v.1.40.0), igraph(v.2.0.3), mime(v.0.12), lifecycle(v.1.0.4), pkgconfig(v.2.0.3), R6(v.2.5.1), fastmap(v.1.2.0), GenomeInfoDbData(v.1.2.12), fitdistrplus(v.1.1-11), future(v.1.33.2), shiny(v.1.8.1.1), digest(v.0.6.36), colorspace(v.2.1-0), tensor(v.1.5), irlba(v.2.3.5.1), progressr(v.0.14.0), fansi(v.1.0.6), spatstat.sparse(v.3.1-0), polyclip(v.1.10-6), httr(v.1.4.7), abind(v.1.4-5), compiler(v.4.4.1), withr(v.3.0.0), BiocParallel(v.1.38.0), MASS(v.7.3-61), DelayedArray(v.0.30.1), rjson(v.0.2.21), tools(v.4.4.1), lmtest(v.0.9-40), httpuv(v.1.6.15), future.apply(v.1.11.2), goftest(v.1.2-3), glue(v.1.7.0), restfulr(v.0.0.15), nlme(v.3.1-165), rhdf5filters(v.1.16.0), promises(v.1.3.0), Rtsne(v.0.17), cluster(v.2.1.6), reshape2(v.1.4.4), generics(v.0.1.3), spatstat.data(v.3.1-2), tidyr(v.1.3.1), sp(v.2.1-4), utf8(v.1.2.4), spatstat.geom(v.3.2-9), RcppAnnoy(v.0.0.22), ggrepel(v.0.9.5), RANN(v.2.6.1), pillar(v.1.9.0), later(v.1.3.2), splines(v.4.4.1), lattice(v.0.22-5), deldir(v.2.0-4), survival(v.3.7-0), tidyselect(v.1.2.1), miniUI(v.0.1.1.1), pbapply(v.1.7-2), knitr(v.1.47), scattermore(v.1.2), xfun(v.0.45), stringi(v.1.8.4), UCSC.utils(v.1.0.0), lazyeval(v.0.2.2), yaml(v.2.3.8), evaluate(v.0.24.0), codetools(v.0.2-20), tibble(v.3.2.1), cli(v.3.6.3), uwot(v.0.2.2), xtable(v.1.8-4), reticulate(v.1.38.0), munsell(v.0.5.1), jquerylib(v.0.1.4), spatstat.random(v.3.2-3), globals(v.0.16.3), png(v.0.1-8), XML(v.3.99-0.17), parallel(v.4.4.1), bitops(v.1.0-7), listenv(v.0.9.1), viridisLite(v.0.4.2), scales(v.1.3.0), ggridges(v.0.5.6), leiden(v.0.4.3.1), purrr(v.1.0.2), crayon(v.1.5.3), rlang(v.1.1.4) and cowplot(v.1.1.3)