CosMx Pancreas Analysis Tutorial

1 Introduction

This tutorial contains the analysis of CosMx spatial transcriptomics data from pancreatic tissue. The analysis includes data processing, visualization, and differential expression analysis.

2 Install libraries

# Install required packages if not already installed
# Remove the hastags if you need to download these

# install.packages("Seurat")
# install.packages("tidyverse")
# install.packages("future")
# install.packages("BiocManager")
# install.packages("openxlsx")
# BiocManager::install("EnhancedVolcano")

3 Load libraries

library(Seurat)
library(tidyverse)
library(openxlsx)
library(ggplot2)
library(future)
library(EnhancedVolcano)
# Sets working environment to hold larger data objects
options(future.globals.maxSize = 8000 * 1024^2)

4 Plot colors

# Define color palette using glasbey and polychrome colors (uses Seurat package)
colors_20 <- DiscretePalette(20, palette = "glasbey", shuffle = FALSE)
colors_10 <- DiscretePalette(10, palette = "polychrome", shuffle = FALSE)

5 Read in data- change readRDS to the path to where your data object is

# This is a Seurat object
## Change path!
so_panc <- readRDS("~/Documents/GCB_Academy_2025/spatial_biology/cosmxsinglecell/so_panc_cosmxWTA_sct.rds") 
Idents(so_panc) <- "types"

6 Basic Visualization

6.0.1 Single FOV (field of view)

# Zoomed to single FOV: CellOverlay_F053.jpg
# Visualize cell types in zoom 1
ImageDimPlot(object = so_panc, fov = "zoom1", group.by = "types", axes = TRUE, cols = "glasbey")

6.0.2 Multiple FOVs

# Zoomed to 6 FOVs: CellOverlay_F051.jpg - CellOverlay_F056.jpg
# Visualize cell types in zoom 4
ImageDimPlot(object = so_panc, fov = "zoom4", group.by = "types", axes = TRUE, cols = "glasbey")

7 Gene Expression Visualizations

7.0.1 Assign marker genes to variables

molsAll <- c("GCG", "EPCAM", "PRSS1", "C1QA",
          "ACTA2", "CD3E", "INS", "CD163", "CD86")
molsDuct <- c("KRT19", "SOX9", "CFTR", "MUC1")
molsBeta <- c("INS", "IAPP", "DLK1", "PDX1")
molsMacs <- c("CD68", "CD163", "CD86", "CSF1R", "MARCO")
molsPSC <- c("FAP", "ACTA2", "COL1A1", "PDGFRA")

7.0.2 Ductal Markers

# Plot ductal markers (molecules) in single zoom segment
ImageDimPlot(so_panc, fov = "zoom1", cols = "polychrome", alpha = 0.25, molecules = molsDuct, 
             mols.size = 0.77, nmols = 20000, border.color = "black", coord.fixed = FALSE)

7.0.3 Beta Cell Markers

# Plot Beta cell markers
ImageDimPlot(so_panc, fov = "zoom1", cols = "polychrome", alpha = 0.25, molecules = molsBeta, 
             mols.size = 0.77, nmols = 20000, border.color = "black", coord.fixed = FALSE)

7.0.4 Macrophage Markers

# Plot macrophage markers
ImageDimPlot(so_panc, fov = "zoom1", cols = "polychrome", alpha = 0.25, molecules = molsMacs, 
             mols.size = 0.77, nmols = 20000, border.color = "black", coord.fixed = FALSE)

7.0.5 PSC (pancreatic stellate cell) Markers

# Plot PSC markers- these are the fibroblasts of the pancreas
ImageDimPlot(so_panc, fov = "zoom1", cols = "polychrome", alpha = 0.25, molecules = molsPSC, 
             mols.size = 0.77, nmols = 20000, border.color = "black", coord.fixed = FALSE)

7.0.6 All Cells Marker Genes

# Plot a group of marker genes for all cell types
ImageDimPlot(so_panc, fov = "zoom1", cols = "polychrome", alpha = 0.25, molecules = molsAll, 
             mols.size = 0.35, nmols = 20000, border.color = "black", coord.fixed = FALSE)

7.1 Cell Type Specific Genes in Specific Cells

# Plot cell type specific genes in specific cells
ImageDimPlot(so_panc, fov = "zoom1", cells = WhichCells(so_panc, 
             # subset to the cells that you want to visualize
             idents = c("Acinar.1", "Acinar.2", "Ductal","Macrophage")), 
             cols = c("red", "green", "blue", "white"), 
             molecules = c("SLC4A4", "SYCN", "LYZ"),
             alpha = 0.3, mols.size = 0.77, mols.cols = c("cyan", "yellow", "purple"))

8 Single-Cell Analysis: all cells in the dataset (all FOVs)

8.0.1 UMAP Visualization - visualize clusters

# Visualize cluster UMAP
DimPlot(so_panc, group.by = "seurat_clusters") + scale_color_manual(values = colors_20)

8.0.2 Visualize cell types

# Visualize cell type UMAP
DimPlot(so_panc) + scale_color_manual(values = colors_20)

8.1 Gene Expression Plots

8.1.1 Violin plot

# Gene expression visualization
VlnPlot(so_panc, features = c("VIM", "PRSS1"), assay = "SCT", 
        layer = "counts", pt.size = 0.1) + NoLegend()

8.1.2 Feature plot

FeaturePlot(so_panc, features = c("KRT17", "VIM", "CD163", "PRSS1"), 
            max.cutoff = "q95")

• What do you notice about the Acinar cell marker, PRSS1?

9 Differential Expression Analysis

9.0.1 Find Markers between Acinar.1 and Acinar.2:

I already ran the FindMarkers algorithm in Seurat since it took me almost 10mins on my laptop. Read in the differentially expressed genes table. Change the path to where you downloaded the data to.

# Find markers between Acinar.1 and Acinar.2 ## This took me ~8mins on my laptop
# Since it takes so long, I provided the output in an Excel file
# 
# acini_markers <- FindMarkers(so_panc, ident.1 = "Acinar.1",
#                             ident.2 = "Acinar.2", verbose = TRUE, 
#                             assay = "SCT")

acini_markers <- openxlsx::read.xlsx("~/Documents/GCB_Academy_2025/spatial_biology/git_cosmx_analysis/Cosmx_WTA/acini_markers.xlsx", 
                    rowNames = TRUE) ### Change the path to where your excel file is!

9.0.2 Volcano Plot of DEGs

# Volcano plot parameters
fc <- 1
pval <- 0.05

degs <- acini_markers %>%
  filter(abs(avg_log2FC) > fc) %>% 
  filter(p_val_adj <= pval) 

# Create volcano plot
EnhancedVolcano(acini_markers,
                lab = rownames(acini_markers),
                selectLab = rownames(degs),
                subtitle = bquote(italic("Acinar.2 <---> Acinar.1")),
                x = 'avg_log2FC',
                y = 'p_val_adj',
                legendLabels = c("NS", paste0("Log2FC >", fc),
                                paste0("Adj. p-val <", pval), 
                                paste0("Adj. p-val <", pval, " and Log2FC >", fc)),
                subtitleLabSize = 15,
                captionLabSize = 15,
                pCutoff = pval,
                FCcutoff = fc,
                pointSize = 2.5,
                labSize = 4,
                axisLabSize = 14,
                legendLabSize = 14, 
                colAlpha = 0.35)

10 Interactive Exploration

10.0.1 See if/where your favorite genes show up in the pancreas

my_genes <- c("SOX9") # Add up to 5 genes here

ImageDimPlot(so_panc, fov = "zoom4", # you can also change this to "zoom1"
             cols = "polychrome", alpha = 0.35, molecules = my_genes, 
             mols.size = 0.35, nmols = 20000, border.color = "black", coord.fixed = FALSE)