Integrate spatial data
11 November, 2025
Load data and libraries
##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(tidyseurat)
library(grid)
library(cowplot)
library(harmony)
source("../bin/plotting_functions.R")
#########
# PATHS #
#########
input_dir <- "../results/01_QC_st_data/"
result_dir <- "../results/02_integrate_st_data/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }
#############
# LODA DATA #
#############
DATA <- readRDS(paste0(input_dir,"seuratObj_filtered.RDS"))Identify Highly Variable Genes (HVG) across samples
################################
# SPLIT INTO SEPERATE DATASETS #
################################
DATA_nested <- DATA %>%
mutate(batch = orig.ident) %>%
nest(data = -batch) %>%
mutate(data = imap(
data, ~ .x %>%
NormalizeData(., normalization.method = "LogNormalize",
verbose = FALSE) %>%
FindVariableFeatures(selection.method = "vst",
nfeatures = 2000,
verbose = FALSE) )) %>%
mutate(data = setNames(.[["data"]], .$batch))
#########################################
# FIND HIGLY VARIABLE GENES PER DATASET #
#########################################
hvgs_heat <- DATA_nested %>%
.$data %>%
map(., ~ .x@assays$RNA@var.features) %>%
# Flatten and collect all HVGs across samples into one vector 'hvgs_all'
( function(x){unique(unlist(x)) ->> hvgs_all; return(x)} ) %>%
# Find intersection of HVGs across samples
( function(x){Reduce(intersect, x) ->> hvgs; return(x)} ) %>%
# Create a df where each column represents a sample, and each row a gene.
# The value is TRUE/FALSE depending on whether that gene is an HVG in that sample.
imap_dfc(., ~hvgs_all %in% .x, .id=.y) %>%
mutate(rownames = hvgs_all) %>%
column_to_rownames(var = "rownames")
# choose the hvg present in at least two samples:
hig_var <- rownames(hvgs_heat)[rowSums(hvgs_heat)>2]
# remove all VDJ-genes from list of HVG
remove <- str_subset(hig_var, "^IGH|^IGK|^IGL|^TRA|^TRB|^TRD|^TRG")
hig_var <- setdiff(hig_var, remove)
length(hig_var)## [1] 4639Heatmap of HVG in all samples
ph <- pheatmap::pheatmap(
t(hvgs_heat * 1),
fontsize_col = 5,
breaks = c(-0.5, 0.5, 1.5),
legend_breaks = c(0, 1), # what values to show in legend
legend_labels = c("No", "Yes"), # how to label them
cluster_rows = FALSE, silent = T,
color = c("grey90", "grey20"),
)
###############################
# SHOW N'TH COLUMN LABELS #
###############################
# Extract column order
col_order <- ph$tree_col$order
labels_ordered <- ph$tree_col$labels[col_order]
# Thin labels: only 1 in 100
labels_thinned <- ifelse(seq_along(labels_ordered) %% 100 == 0, ph$tree_col$labels, "") # show every 100th label
# Replace column labels in the gtable (this actually affects the plot)
ph[["gtable"]][["grobs"]][[3]][["label"]] <- labels_thinned
# Draw the modified heatmap
#grid.newpage()
grid.draw(ph$gtable)
Showing every 100th column label of the heatmap
Integration
############
# HARMONY #
###########
DATA <- DATA %>%
# SCTransform(verbose = FALSE, ) %>%
NormalizeData(verbose = FALSE) %>%
FindVariableFeatures(selection.method = "vst",
nfeatures = 4000,
verbose = FALSE) %>%
ScaleData(verbose = FALSE, features = hig_var ) %>%
RunPCA(verbose = FALSE, npcs = 50) %>%
RunUMAP(dims = 1:50,
n.components = 2L,
n.neighbors = 10,
min.dist = .1,
spread = .3)
DATA <- DATA %>%
RunHarmony(group.by.vars = "orig.ident",
reduction.use = "pca",
dims.use = 1:50,
assay.use = "RNA") #%>%
DATA <- DATA %>%
RunUMAP(dims = 1:50,
n.neighbors = 10,
min.dist = .1,
spread = 1,
repulsion.strength = 1,
negative.sample.rate = 10,
n.epochs = 100,
reduction = "harmony",
reduction.name = "umapharmony")Alternative graph based UMAP
integrated <- DATA@reductions$harmony@cell.embeddings
ann <- RcppHNSW::hnsw_build(as.matrix(integrated), distance = "cosine")
knn <- RcppHNSW::hnsw_search(as.matrix(integrated) , ann = ann , k = 15)
UU2 <- uwot::umap(X = NULL,
nn_method = knn,
n_components = 2,
ret_extra = c("model","fgraph"),
verbose = T,
min_dist = 0.1,
spread = .3,
repulsion_strength = 1,
negative_sample_rate = 10,
n_epochs = 150,
n_threads = 8)
dimnames(UU2$embedding) <- list(colnames(DATA),paste0("umap_harmony_knn_", 1:2))
DATA@reductions[["umap_harmony_knn"]] <- CreateDimReducObject(embeddings = UU2$embedding,
key = "umap_harmony_knn_")
colnames(DATA@reductions$umap_harmony_knn@cell.embeddings) <- paste0("umap_harmony_knn_", 1:2)res <- c("umapharmony", "umap_harmony_knn")
p <- map(res, ~plot_clusters.fun(DATA, red=.x, cluster="orig.ident", lable=FALSE, txt_size = 7))
plot_grid(ncol = 2,
plotlist = p)
Plot before and after integration
# dev.new(height=6, width=6.6929133858, noRStudioGD = TRUE)
res <- c("PC", "harmony", "UMAP", "umapharmony")
title <- c("PCA raw data", "PCA Harmony integrated", "UMAP raw data", "UMAP Harmony integrated")
p <- map2(res, title,
~plot_clusters.fun(DATA,
cluster="orig.ident", txt_size = 9,
red=.x, lable=FALSE, title=.y))
plot_grid(ncol = 2,
plotlist = p)