Reproducing the figures of the manuscript

Introduction

This markdown document reproduces the figures of the manuscript: Multi-omic identification of perineurial hyperplasia and lipid-associated nerve macrophages in human polyneuropathies by Heming et al, Nature Communications, August 2025.

All relevant data are automatically downloaded from Zenodo.

Instructions on how to restore the environment using renv or Docker can be found in the corresponding GitHub repository. The corresponding website, which contains interactive visualizations of the single nucleus transcriptomics and spatial transcriptomics data, is available at pns-atlas.mzhlab.com. If you have any questions, please contact us at mheming.com.

Load libraries

library(qs)
library(Seurat)
library(tidyverse)
library(miloDE)
library(sessioninfo)
library(EnhancedVolcano)
library(pals)
library(scMisc)
library(liana)

Figure 1

UMAP of all cells (Figure 1B)

options(timeout = 3600)
if (!file.exists("umap_figure.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/umap_figure.qs?download=1",
        "umap_figure.qs"
    )
}
umap_figure <- qread("umap_figure.qs")
DimPlot(
    umap_figure,
    reduction = "umap.scvi.full",
    pt.size = .1,
    alpha = .1,
    cols = umap_figure@misc$cluster_col,
    label = TRUE,
    raster = FALSE
) +
    NoLegend() +
    theme(
        axis.text = element_blank(),
        axis.ticks = element_blank()
    ) +
    xlab("UMAP_1") +
    ylab("UMAP_2")

Figure 2

UMAP immune cells (Figure 2A)

if (!file.exists("ic_figure.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/ic_figure.qs?download=1",
        "ic_figure.qs"
    )
}
ic_figure <- qread("ic_figure.qs")
DimPlot(
    ic_figure,
    reduction = "umap.rpca",
    pt.size = .1,
    alpha = .3,
    cols = ic_figure@misc$ic_cluster_col,
    label = TRUE,
    raster = FALSE
) +
    theme(
        axis.text = element_blank(),
        axis.ticks = element_blank()
    ) +
    NoLegend() +
    xlab("UMAP1") +
    ylab("UMAP2")

Feature plots immune cells (Figure 2B)

FeaturePlot(
    ic_figure,
    features = c("MS4A7", "CX3CR1", "TREM2", "LYVE1", "FOLR2", "TIMD4"),
    reduction = "umap.rpca",
    pt.size = 0.1,
    raster = FALSE,
    coord.fixed = TRUE,
    cols = c("#F0F0F0", "#CB181D"),
    order = TRUE,
    ncol = 3
) &
    theme(
        axis.text = element_blank(),
        axis.ticks = element_blank()
    ) &
    xlab("UMAP1") &
    ylab("UMAP2")

Enrichment of gene ontology terms in markers genes of Macro 18 (Figure 2B)

if (!file.exists("enrichr_macro18.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/enrichr_macro18.qs?download=1",
        "enrichr_macro18.qs"
    )
}
enrichr_macro18 <- qread("enrichr_macro18.qs")
enrichr_macro18 |>
    dplyr::slice_min(order_by = Adjusted.P.value, n = 10, with_ties = FALSE) |>
    tidyr::separate(Overlap, into = c("overlap1", "overlap2")) |> # separate overlap in two columns
    dplyr::mutate(
        Term = gsub(x = Term, pattern = "\\s\\(.+\\)", replacement = "")
    ) |>
    dplyr::mutate(overlap = as.numeric(overlap1) / as.numeric(overlap2)) |> # calculcate overlap
    ggplot(aes(
        y = reorder(Term, -log10(Adjusted.P.value)),
        x = -log10(Adjusted.P.value)
    )) +
    geom_col(fill = scales::hue_pal()(5)[1]) +
    labs(
        x = "-Log10 Adjusted P value",
        y = ""
    ) +
    theme_classic() +
    theme(legend.position = "none")

Dotplot of SAMC markers in IC (Figure 2D)

if (!file.exists("dotplot_data_ic_samc.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/dotplot_data_ic_samc.qs?download=1",
        "dotplot_data_ic_samc.qs"
    )
}

# Seurat DotPlot function (only ggplot part to avoid calculations)
DotPlotModified <- function(
    data.plot,
    cols = c("lightgrey", "blue"),
    dot.scale = 6,
    scale.by = "radius",
    scale.min = NA,
    scale.max = NA) {
    scale.func <- switch(
        EXPR = scale.by,
        "size" = scale_size,
        "radius" = scale_radius,
        stop("'scale.by' must be either 'size' or 'radius'")
    )
    color.by <- "avg.exp.scaled"
    plot <- ggplot(
        data = data.plot,
        mapping = aes_string(x = "features.plot", y = "id")
    ) +
        geom_point(mapping = aes_string(size = "pct.exp", color = color.by)) +
        scale.func(range = c(0, dot.scale), limits = c(scale.min, scale.max)) +
        theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
        guides(size = guide_legend(title = "Percent Expressed")) +
        cowplot::theme_cowplot()
    plot <- plot + guides(color = guide_colorbar(title = "Average Expression"))
    return(plot)
}

dotplot_data_ic_samc <- qread("dotplot_data_ic_samc.qs")

DotPlotModified(
    data.plot = dotplot_data_ic_samc,
    scale.by = "size"
) +
    viridis::scale_color_viridis(option = "viridis") +
    theme(
        axis.text.x = element_text(
            angle = 90,
            vjust = 0.5,
            hjust = 1,
            face = "italic"
        ),
    ) +
    xlab("") +
    ylab("")

Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.

Warning: Removed 17 rows containing missing values or values outside the scale range
(`geom_point()`).

Feature plots of SAMC markers in IC (Figure 2E)

ic_figure <- qread("ic_figure.qs")
plot_predicted_samc <- function(seu_obj) {
    # Create a data frame with alpha values, color values, and UMAP coordinates
    umap_coords <- Embeddings(seu_obj, "umap.rpca")
    alpha_values <- ifelse(seu_obj$stroke_label == "SAMC", 0.5, 0.001)
    color_values <- ifelse(seu_obj$stroke_label == "SAMC", "SAMC", "Other")
    color_map <- setNames(c("red", "grey"), c("SAMC", "Other"))
    plot_df <- data.frame(
        cell = rownames(umap_coords),
        UMAP1 = umap_coords[, 1],
        UMAP2 = umap_coords[, 2],
        alpha = alpha_values,
        color = color_values
    )

    # Create custom ggplot
    samc_predicted <- ggplot(
        plot_df,
        aes(x = UMAP1, y = UMAP2, alpha = alpha, color = color)
    ) +
        geom_point(size = 0.1) +
        scale_color_manual(values = color_map) + # Simplified color scale
        theme_minimal() +
        theme(
            axis.text = element_blank(),
            axis.ticks = element_blank(),
            panel.grid = element_blank(),
            panel.border = element_rect(color = "black", fill = NA, size = 1.0),
            panel.background = element_rect(fill = "white", color = NA),
            plot.background = element_rect(fill = "white", color = NA),
            legend.position = "none"
        ) +
        ggtitle("predicted SAMC")
    return(samc_predicted)
}

plot_predicted_samc(ic_figure)

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.

Figure 3

Propeller abundance PNP vs CTRL (Figure 3A)

if (!file.exists("propeller_PNP_CTRL.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/propeller_PNP_CTRL.qs?download=1",
        "propeller_PNP_CTRL.qs"
    )
}
if (!file.exists("propeller_PNP_CTRL_ic.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/propeller_PNP_CTRL_ic.qs?download=1",
        "propeller_PNP_CTRL_ic.qs"
    )
}
propeller_PNP_CTRL <- qread("propeller_PNP_CTRL.qs")
propeller_PNP_CTRL_ic <- qread("propeller_PNP_CTRL_ic.qs")

# function to plot custom dotplot
dotplotPropeller <- function(data, color, title) {
    ggplot(
        data,
        aes(x = log2ratio, y = fct_reorder(cluster, log2ratio), color = cluster)
    ) +
        geom_point(size = 5) +
        theme_classic() +
        geom_vline(
            xintercept = 0,
            color = "red",
            linetype = "solid"
        ) +
        scale_color_manual(values = color) +
        xlab("Log2 fold change") +
        ylab(NULL) +
        theme(legend.position = "none") +
        ggtitle(title)
}
dotplotPropeller(
    propeller_PNP_CTRL,
    color = umap_figure@misc$cluster_col,
    title = "PNP vs CTRL main clusters"
)

dotplotPropeller(
    propeller_PNP_CTRL_ic,
    color = ic_figure@misc$ic_cluster_col,
    title = "PNP vs CTRL immune cells"
)

Covarying Neighborhood Analysis (Figure 3B)

if (!file.exists("cna_pnp_gratio_figure.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/cna_pnp_gratio_figure.qs?download=1",
        "cna_pnp_gratio_figure.qs"
    )
}
if (!file.exists("cna_incat_figure.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/cna_incat_figure.qs?download=1",
        "cna_incat_figure.qs"
    )
}
cna_pnp_gratio_figure <- qread("cna_pnp_gratio_figure.qs")
cna_incat_figure <- qread("cna_incat_figure.qs")

# function to plot CNA feature plot
fplotCNA <- function(object, feature, title) {
    FeaturePlot(
        object,
        reduction = "umap.scvi.full",
        features = feature,
        pt.size = 0.1,
        order = FALSE,
        coord.fixed = TRUE,
        raster = FALSE,
        alpha = 0.2
    ) +
        scale_colour_gradient2(
            low = "#2166AC",
            mid = "white",
            high = "#B2182B",
            midpoint = 0,
        ) +
        theme(
            axis.text = element_blank(),
            axis.ticks = element_blank(),
        ) +
        labs(
            title = title,
            color = "Correlation",
            x = "UMAP1",
            y = "UMAP2"
        )
}

fplotCNA(cna_pnp_gratio_figure, "cna_ncorrs_pnp", title = "PNP")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

fplotCNA(cna_incat_figure, "cna_ncorrs", title = "INCAT")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

fplotCNA(cna_pnp_gratio_figure, "cna_ncorrs_gratio", title = "g-ratio")

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Number of DEGs (Figure 3C)

if (!file.exists("pnp_ctrl_pseudobulk_de.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/pnp_ctrl_pseudobulk_de.qs?download=1",
        "pnp_ctrl_pseudobulk_de.qs"
    )
}
pnp_ctrl_pseudobulk <- qread("pnp_ctrl_pseudobulk_de.qs")
pnp_ctrl_pseudobulk |>
    mutate(cluster = fct_reorder(cluster, n)) |>
    ggplot(aes(x = cluster, y = n, fill = cluster)) +
    geom_col() +
    coord_flip() +
    scale_fill_manual(values = umap_figure@misc$cluster_col) +
    theme_classic() +
    theme(legend.position = "none") +
    labs(
        x = "",
        y = "",
        title = "PNP vs CTRL"
    )

MiloDE DEGs PNP vs CTRL (Figure 3D)

if (!file.exists("milo_figure.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/milo_figure.qs?download=1",
        "milo_figure.qs"
    )
}

milo_figure <- qread("milo_figure.qs")
plotMiloDE <- function(condition, title) {
    plot <-
        plot_milo_by_single_metric(
            milo_figure$obj,
            milo_figure$stat[[condition]],
            colour_by = "n_DE_genes",
            layout = "UMAP.SCVI.FULL",
            size_range = c(0.5, 5),
            edge_width = c(0.1, 1.0),
            edge_weight.thres = 10
        ) +
        viridis::scale_fill_viridis(name = "# DE genes", option = "inferno") +
        ggtitle(title)
    print(plot)
}

plotMiloDE("pnp", "PNP vs CTRL")

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Pseudobulk DE PNP vs CTRL (Figure 3E)

if (!file.exists("de_pseudo_pnp_ctrl.qs")) {
    download.file(
        "https://zenodo.org/records/15108216/files/de_pseudo_pnp_ctrl.qs?download=1",
        "de_pseudo_pnp_ctrl.qs"
    )
}
de_pseudo_pnp_ctrl <- qread("de_pseudo_pnp_ctrl.qs")

# volcano plot function
volcanoPlot <- function(
    cluster,
    input,
    FCcutoff = 2,
    selectLab = NULL,
    drawConnectors = TRUE,
    condition1,
    condition2) {
    input <- input[[cluster]]
    volcano <- EnhancedVolcano::EnhancedVolcano(
        data.frame(input),
        lab = paste0("italic('", input[["gene"]], "')"),
        x = "avg_logFC",
        y = "p_val_adj",
        xlim = c(min(input[["avg_logFC"]], max(input[["avg_logFC"]]))),
        ylim = c(0, max(-log10(input[["p_val_adj"]]))),
        pCutoff = 0.1,
        FCcutoff = FCcutoff,
        axisLabSize = 15,
        pointSize = 2,
        labSize = 5,
        subtitle = NULL,
        caption = NULL,
        border = "full",
        gridlines.major = FALSE,
        gridlines.minor = FALSE,
        drawConnectors = drawConnectors,
        lengthConnectors = unit(0.0001, "npc"),
        title = paste(condition1, "vs", condition2, "in ", cluster),
        boxedLabels = TRUE,
        selectLab = selectLab,
        xlab = bquote(~ Log[2] ~ "fold change"),
        ylab = bquote(~ -Log[10] ~ "adjusted p-value"),
        parseLabels = TRUE,
        legendLabels = c(
            "NS",
            "logFC",
            "p-val",
            "p-val + logFC"
        ),
        legendPosition = "right",
    )
}

# define clusters of interest
cluster_de <- c("mySC", "nmSC", "repairSC", "PC2")

# define genes of interests
lab_pnp_ctrl <- list(
    "mySC" = paste0(
        "italic('",
        c("DCN", "TNXB", "COL1A1", "COL15A1", "CD53", "IL4R", "CD74"),
        "')"
    ),
    "nmSC" = paste0(
        "italic('",
        c("IL10RA", "IL13RA1", "CSF2RA", "TGFBI"),
        "')"
    ),
    "repairSC" = paste0("italic('", c("GALR1", "TMEM47"), "')"),
    "PC2" = paste0(
        "italic('",
        c("MFAP5", "NLGN4Y", "PCDH11Y", "IFIT3", "OASL", "MX1"),
        "')"
    )
)

# plot volcano plots of selected clusters
purrr::walk(
    cluster_de,
    function(cluster) {
        print(volcanoPlot(
            cluster = cluster,
            input = de_pseudo_pnp_ctrl,
            FCcutoff = 2,
            condition1 = "PNP",
            condition2 = "CTRL",
            selectLab = lab_pnp_ctrl[[cluster]]
        ))
    }
)