Figure 1 (B)
# Load the library
library(aplotExtra)
library(tidyverse)
library(ggfun)
library(readxl)
# Define functions and parameters
common_theme <- theme(
legend.position = "none",
axis.text.x = element_text(angle = 20, hjust = 1, size = 12, colour = "black")
)
params_tbl <- tribble(
~type, ~d_id, ~col, ~low, ~high,
"bar", "d2", 2, "#FFFFE5", "#cc4c02",
"point", "d1", 3, "#FFFFE5", "#cc4c02",
"bar", "d2", 4, "#D0E9C6", "#4CAF50",
"bar", "d2", 6, "#D0E9C6", "#4CAF50",
"point", "d1", 7, "#D0E9C6", "#4CAF50",
"point", "d1", 8, "#BBDEFB", "#2196F3",
"point", "d1", 9, "#BBDEFB", "#2196F3",
"point", "d1", 10, "#BBDEFB", "#2196F3",
"bar", "d2", 8, "#CAC0E1", "#715EA9",
"bar", "d2", 10, "#CAC0E1", "#715EA9",
"point", "d1", 14, "#CAC0E1", "#715EA9",
"point", "d1", 15, "#CAC0E1", "#715EA9"
)
custom_point <- function(d, col, low, high) {
funky_point(d, col) +
scale_size_continuous(limits = c(0, 1), range = c(1, 5)) +
scale_fill_gradient(low = low, high = high) +
common_theme
}
custom_bar <- function(d, col, low, high) {
p <- funky_bar(d, col)
p$layers <- p$layers[-length(p$layers)]
p +
geom_vline(xintercept = 1.2) +
scale_fill_gradient(low = low, high = high) +
common_theme
}
plot_heatmap <- function(d1, d2, params = params_tbl) {
f1 <- funky_text(d1)
plots_list <- pmap(params, function(type, d_id, col, low, high) {
df <- if(d_id == "d1") d1 else d2
if(type == "bar") {
custom_bar(df, col, low, high)
} else {
custom_point(df, col, low, high)
}
})
ordered <- c(list(f1), plots_list)
do.call(funky_heatmap, c(ordered, list(
options = list(
theme_stamp(),
ggfun::theme_blinds(
colour = c('white','white','grey90','grey90'),
axis = 'y',
axis.line.x = ggplot2::element_line()
)
)
)))
}
# Load the data and plot figure
d1_file <- "heatmap1030.xlsx"
d2_file <- "heatmap1031.xlsx"
stopifnot(file.exists(d1_file), file.exists(d2_file))
d1 <- read_excel(d1_file) |> as.data.frame() |> tibble::column_to_rownames("Sample_name")
d2 <- read_excel(d2_file) |> as.data.frame()
plot_heatmap(d1, d2)
Figure 1 (C)
# Load the library
library(Seurat)
library(ggprism)
#set the color
custom_colors <- c(
"Brain1 NV" = "#ED9F9B",
"Brain1 SF" = "#E5A79A",
"Brain2 NV" = "#F1D0C6",
"Brain2 SF" = "#6E8FB2",
"Brain3 NV" = "#8AB1D2",
"Brain3 SF" = "#A6CDE4",
"Brain4 NV" = "#D0E0EF",
"Brain4 SF" = "#74B69F",
"Brain5 NV" = "#80BA8A",
"Brain5 SF" = "#9CD1C8",
"Lung NV" = "#F1DFA4",
"Lung SF" = "#CC88B0"
)
# Load the data
load("merged_10Brain_Lung.RData")
data <- merged_10Brain_Lung@meta.data[, c("nCount_RNA", "nFeature_RNA", "orig.ident")]
# Show the figure
p <- ggplot(data, aes(x = nCount_RNA, y = nFeature_RNA, color = orig.ident)) +
stat_smooth(method = "loess", se = FALSE, size = 1.2) +
labs(x = "number of UMI", y = "number of genes") +
scale_color_manual(values = custom_colors) +
theme_minimal() +
theme_prism(base_size = 15) +
theme(
plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
axis.title = element_text(size = 18, face = "bold"),
axis.text = element_text(size = 15, face = "bold"),
legend.text = element_text(size = 12, face = "bold"),
plot.margin = margin(25, 25, 25, 25)
) +
guides(color = guide_legend(
title = "Platform",
title.theme = element_text(size = 13, face = "bold"),
label.theme = element_text(size = 12, face = "bold")
))
p
Figure 1 (D)
# Load the library
library(eulerr)
library(grid)
# Define functions and parameters
format_labels <- function(counts, digits = 4) {
total <- sum(counts)
if (total == 0) {
out <- rep("0\n(0%)", length(counts))
names(out) <- names(counts)
return(out)
}
out <- paste0(counts, "\n(", signif(counts / total * 100, digits), "%)")
names(out) <- names(counts)
out
}
extract_barcodes <- function(seurat_obj, ident) {
md <- seurat_obj@meta.data
stopifnot("orig.ident" %in% colnames(md))
cells <- rownames(md)[md$orig.ident == ident]
sub(paste0("^", ident, "_"), "", cells)
}
analyze_barcode_overlap <- function(seurat_obj, sample1, sample2, percentage_digits = 4) {
b1 <- unique(extract_barcodes(seurat_obj, sample1))
b2 <- unique(extract_barcodes(seurat_obj, sample2))
shared_n <- length(intersect(b1, b2))
counts <- c(
setNames(length(setdiff(b1, b2)), sample1),
setNames(length(setdiff(b2, b1)), sample2),
setNames(shared_n, paste(sample1, sample2, sep = "&"))
)
labels <- format_labels(counts, percentage_digits)
list(data = counts, labels = labels)
}
plot_barcode_euler <- function(overlap_res, title = "Sample") {
venn <- euler(overlap_res$data)
venn_labels <- overlap_res$labels
grid.newpage()
p <- plot(
venn,
fills = list(fill = c("#F9DF91", "#A8CFE8", "#E6D9AE"), alpha = 0.7),
labels = FALSE,
edges = NULL,
quantities = FALSE,
main = list(label = title, cex = 4, font = 2),
output = "grob"
)
grid.draw(p)
grid.text(venn_labels[1], x = 0.12, y = 0.5, gp = gpar(fontsize = 20, fontface = "bold"))
grid.text(venn_labels[2], x = 0.88, y = 0.5, gp = gpar(fontsize = 20, fontface = "bold"))
grid.text(venn_labels[3], x = 0.5, y = 0.5, gp = gpar(fontsize = 40, fontface = "bold"))
}
# Load data and show the figure
load("merged_all.RData")
res <- analyze_barcode_overlap(merged_R_9_24h_2, "Lung NV", "Lung SF")
plot_barcode_euler(res, title = "Lung")
Figure 1 (E)
# Define functions and parameters
plot_params <- list(
gene_col = "nFeature_RNA",
umi_col = "nCount_RNA",
bar_colors = c("#F9DF91", "#E6D9AE", "#A8CFE8"),
show_overall_median = TRUE
)
get_barcodes <- function(meta, ident) {
cells <- rownames(meta)[meta$orig.ident == ident]
sub(paste0(ident, "_"), "", cells)
}
prepare_data <- function(seurat_obj, sample1, sample2, gene_col, umi_col) {
meta <- seurat_obj@meta.data
barcodes1 <- get_barcodes(meta, sample1)
barcodes2 <- get_barcodes(meta, sample2)
shared <- intersect(barcodes1, barcodes2)
unique1 <- setdiff(barcodes1, shared)
unique2 <- setdiff(barcodes2, shared)
cells_sample1 <- paste0(sample1, "_", unique1)
cells_sample2 <- paste0(sample2, "_", unique2)
cells_shared <- c(paste0(sample1, "_", shared), paste0(sample2, "_", shared))
df <- data.frame(
Cell_Group = factor(
c(rep(paste0(sample1, " unique"), length(cells_sample1)),
rep("Shared", length(cells_shared)),
rep(paste0(sample2, " unique"), length(cells_sample2))),
levels = c(paste0(sample1, " unique"), "Shared", paste0(sample2, " unique"))
),
Gene_Count = c(meta[cells_sample1, gene_col],
meta[cells_shared, gene_col],
meta[cells_sample2, gene_col]),
UMI_Count = c(meta[cells_sample1, umi_col],
meta[cells_shared, umi_col],
meta[cells_sample2, umi_col])
)
gene_medians <- tapply(df$Gene_Count, df$Cell_Group, median)
gene_se <- tapply(df$Gene_Count, df$Cell_Group, function(x) mad(x) / sqrt(length(x)))
overall_median <- median(df$Gene_Count)
return(list(df = df, gene_medians = gene_medians,
gene_se = gene_se, overall_median = overall_median))
}
plot_count_barplot <- function(gene_medians, gene_se, overall_median,
bar_colors, show_overall_median) {
par(mar = c(5, 8, 5, 5))
bp <- barplot(
gene_medians,
beside = TRUE,
col = bar_colors,
ylab = "Median gene count per cell",
cex.main = 1.6, cex.axis = 1.2, cex.lab = 1.2,
las = 2, font.lab = 2, font.axis = 2, font.main = 2,
names.arg = FALSE
)
text(
x = bp,
y = par("usr")[3] - 0.05 * diff(par("usr")[3:4]),
labels = names(gene_medians),
srt = 25, adj = 1, xpd = TRUE,
cex = 1.1, font = 2
)
if (show_overall_median) {
abline(h = overall_median, lty = 2, col = "red", lwd = 2)
}
arrows(
x0 = bp, y0 = gene_medians - gene_se,
x1 = bp, y1 = gene_medians + gene_se,
angle = 90, code = 3, length = 0.1, lwd = 1.3
)
}
# Show the figure
res <- prepare_data(
seurat_obj = merged_R_9_24h_2,
sample1 = "Lung NV",
sample2 = "Lung SF",
gene_col = plot_params$gene_col,
umi_col = plot_params$umi_col
)
plot_count_barplot(
gene_medians = res$gene_medians,
gene_se = res$gene_se,
overall_median = res$overall_median,
bar_colors = plot_params$bar_colors,
show_overall_median = plot_params$show_overall_median
)
# package version
library(sessioninfo)
pi <- package_info(pkgs = "attached")
pi_df <- as.data.frame(pi)[, c("package", "loadedversion", "source")]
pi_df
## package loadedversion source
## aplotExtra aplotExtra 0.0.4 CRAN (R 4.3.3)
## dplyr dplyr 1.1.4 CRAN (R 4.3.3)
## eulerr eulerr 7.0.2 CRAN (R 4.3.3)
## forcats forcats 1.0.0 CRAN (R 4.3.1)
## ggfun ggfun 0.1.8 CRAN (R 4.3.3)
## ggplot2 ggplot2 3.5.2 CRAN (R 4.3.3)
## ggprism ggprism 1.0.6 CRAN (R 4.3.3)
## ggstar ggstar 1.0.4 CRAN (R 4.3.1)
## lubridate lubridate 1.9.4 CRAN (R 4.3.3)
## purrr purrr 1.0.4 CRAN (R 4.3.3)
## readr readr 2.1.5 CRAN (R 4.3.1)
## readxl readxl 1.4.5 CRAN (R 4.3.3)
## sessioninfo sessioninfo 1.2.3 CRAN (R 4.3.3)
## Seurat Seurat 5.1.0 CRAN (R 4.3.3)
## SeuratObject SeuratObject 5.1.0 CRAN (R 4.3.3)
## sp sp 2.2-0 CRAN (R 4.3.3)
## stringr stringr 1.5.1 CRAN (R 4.3.1)
## tibble tibble 3.3.0 CRAN (R 4.3.3)
## tidyr tidyr 1.3.1 CRAN (R 4.3.1)
## tidyverse tidyverse 2.0.0 CRAN (R 4.3.3)