Figure 4 (A)
# Load the library
library(eulerr)
library(grid)
# Define functions
compute_venn <- function(nv_genes, sf_genes, sample1 = "Brain NV", sample2 = "Brain SF", digits = 4) {
set1_only <- setdiff(nv_genes, sf_genes)
set2_only <- setdiff(sf_genes, nv_genes)
shared <- intersect(nv_genes, sf_genes)
counts <- c(length(set1_only), length(set2_only), length(shared))
names(counts) <- c(sample1, sample2, paste(sample1, sample2, sep = "&"))
total <- sum(counts)
labels <- paste0(counts, "\n(", signif(counts / total * 100, digits), "%)")
names(labels) <- names(counts)
list(counts = counts, labels = labels)
}
plot_marker_euler <- function(venn_res, title = "Sample", fills = c("#F9DF91", "#E6D9AE", "#A8CFE8")) {
venn_obj <- euler(venn_res$counts)
grid.newpage()
p <- plot(
venn_obj,
fills = list(fill = fills, alpha = 0.7),
labels = FALSE,
edges = NULL,
quantities = FALSE,
main = list(label = title, cex = 2.1, font = 2),
output = "grob"
)
grid.draw(p)
grid.text(venn_res$labels[1], x = 0.14, y = 0.5, gp = gpar(fontsize = 16, fontface = "bold"))
grid.text(venn_res$labels[2], x = 0.86, y = 0.5, gp = gpar(fontsize = 16, fontface = "bold"))
grid.text(venn_res$labels[3], x = 0.5, y = 0.5, gp = gpar(fontsize = 20, fontface = "bold"))
}
# Load data and show the figure
load("marker_10bnvsf.RData")
celltype <- "Astrocytes"
nv_genes <- markers_10bnv$gene[markers_10bnv$cluster == celltype]
sf_genes <- markers_10bsf$gene[markers_10bsf$cluster == celltype]
venn_res <- compute_venn(nv_genes, sf_genes)
plot_marker_euler(venn_res, title = celltype)
Figure 4 (B)
# Load the library
library(Seurat)
library(ggplot2)
library(org.Mm.eg.db)
library(dplyr)
library(clusterProfiler)
# Define functions
prepare_GO <- function(marker_data, p_cutoff = 0.05, q_cutoff = 0.2) {
markers <- marker_data %>%
group_by(cluster) %>%
filter(p_val_adj < 0.001) %>%
ungroup()
gid <- bitr(unique(markers$gene), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
markers <- full_join(markers, gid, by = c("gene" = "SYMBOL"))
compareCluster(
ENTREZID ~ cluster,
data = markers,
fun = "enrichGO",
OrgDb = org.Mm.eg.db,
ont = "ALL",
keyType = "ENTREZID",
pvalueCutoff = p_cutoff,
qvalueCutoff = q_cutoff
)
}
clean_go_description <- function(df) {
df$Description <- gsub(" - Mus musculus \\(house mouse\\)", "", df$Description)
df
}
combine_results <- function(go_obj_list, sample_names) {
combined <- lapply(seq_along(go_obj_list), function(i) {
df <- clean_go_description(go_obj_list[[i]]@compareClusterResult)
df$sample <- sample_names[i]
df
}) %>% bind_rows()
combined_go <- go_obj_list[[1]]
combined_go@compareClusterResult <- combined
combined_go
}
# Load the data and combine
b1_5nv_GO <- prepare_GO(markers_10bnv, p_cutoff = 0.03)
b1_5sf_GO <- prepare_GO(markers_10bsf, p_cutoff = 0.03)
b1_5_combined_GO <- combine_results(
go_obj_list = list(b1_5nv_GO, b1_5sf_GO),
sample_names = c("Brain1-5 NV", "Brain1-5 SF")
)
b1_5_combined_GO@compareClusterResult$cluster <- factor(
b1_5_combined_GO@compareClusterResult$cluster,
levels = unique(b1_5_combined_GO@compareClusterResult$cluster)
)
# Load the data and show the figure
p <- dotplot(
b1_5_combined_GO,
x = "cluster",
color = "p.adjust",
font.size = 18,
showCategory = 5,
label_format = 40
) +
theme(
axis.text.x = element_text(angle = 30, hjust = 1, size = 16, face = "bold"),
axis.text.y = element_text(size = 16, face = "bold"),
axis.title = element_blank(),
plot.title = element_text(size = 32, face = "bold"),
legend.title = element_text(size = 20, face = "bold"),
legend.text = element_text(size = 20, face = "bold"),
strip.text = element_text(size = 22, face = "bold")
) +
facet_grid(. ~ sample)
p
Figure 4 (D)
import warnings
warnings.filterwarnings('ignore')
import stereo as st
import pickle
import matplotlib.pyplot as plt
# Define functions
class SafeUnpickler(pickle.Unpickler):
def find_class(self, module, name):
if module == 'stereo.core.result' and name == 'Result':
class DummyResult(dict):
def __setitem__(self, key, value):
super().__setitem__(key, value)
return DummyResult
return super().find_class(module, name)
def load_stereo_data(file_path):
with open(file_path, "rb") as f:
return SafeUnpickler(f).load()
def subset_cells(data, prefix):
cells = [c for c in data.cells.cell_name if c.startswith(prefix)]
sub_data = data.tl.filter_cells(cell_list=cells, inplace=False)
sub_module = data.tl.result["spatial_hotspot"].module_scores.loc[cells, :]
sub_data.tl.result["spatial_hotspot"].module_scores = sub_module
return sub_data
def plot_module(data_obj, module_id, vmin=-1.2, vmax=2.6, cbar_ticks=None, title=None):
if cbar_ticks is None:
cbar_ticks = [vmin, (vmin+vmax)/2, vmax]
scores = data_obj.tl.result["spatial_hotspot"].module_scores[module_id]
coords = data_obj.cells_matrix["spatial"]
plt.figure(figsize=(6, 5))
sc = plt.scatter(coords[:, 0], coords[:, 1],
c=scores, cmap="RdYlBu_r", s=5,
vmin=vmin, vmax=vmax)
ax = plt.gca()
ax.invert_yaxis()
ax.set_aspect("equal")
ax.set_xticks([]); ax.set_yticks([])
for spine in ax.spines.values():
spine.set_visible(False)
cbar = plt.colorbar(sc, shrink=0.5, aspect=17)
cbar.set_ticks(cbar_ticks)
plt.axis("off")
plt.show()
# Load the data and plot figure
data = load_stereo_data("Brain_NVSF.pkl")
plot_module(subset_cells(data, "Brain1 NV"), module_id=4,
vmin=-1.2, vmax=2.6,
cbar_ticks=[-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5])
# Show the figure
plot_module(subset_cells(data, "Brain1 SF"), module_id=4,
vmin=-1.2, vmax=2.6,
cbar_ticks=[-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5])
