Home

“I see you haven’t changed at all, Captain, in these four years,” Nikolai Vsevolodovich said, as if somewhat more kindly. “It must be true that the whole second half of a man’s life is most often made up only of habits accumulated during the first half.”

Fyodor Dostoevsky, Demons

The purpose of this markdown is to take Seurat’s FeaturePlot function and give it a rainbow color scheme, of the type that we are used to in the cytometry world.

First, let’s get some data and do the standard Seurat workup.

library(tidyverse)
library(Seurat)
library(SeuratData)

We will next pull the PBMC 3k dataset from the SeuratData package.

pbmc <- SeuratData::LoadData(ds = "pbmc3k")
## Warning: Assay RNA changing from Assay to Assay
## Warning: Assay RNA changing from Assay to Assay5

We will next do the standard Seurat workup.

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
pbmc <- ScaleData(pbmc, features = rownames(pbmc))
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- RunUMAP(pbmc, dims = 1:10)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

Now we get to where we want to be. Have a look at the built-in feature plot here.

genes <- c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A")
FeaturePlot(pbmc, features = genes)

This is great and all, but there is a specific color scheme that people who started with CyTOF prefer. I’m going to call this the Cytobank color scheme. And that is more of a rainbow scheme that starts with blue, then runs through cyan, green, yellow, orange, and then red.

The way you get it has some interesting notation. To give credit where credit is due, the solution with the “&” comes from the following discussion thread here.

cytobank_colors <- colorRampPalette(c("blue", "cyan", "green", "yellow", "orange", "red"))

FeaturePlot(pbmc, features = genes) & 
  scale_colour_gradientn(colours = cytobank_colors(100))

You can use whatever scheme you like, but to my CyTOF-tinted eyes, this works for me a bit better. You can see that in some instances, like LYZ, it helps to distinguish between low and high expression. In other instances like CD3E, there is not so much variation on the T cells, so it doesn’t seem to matter as much.

Now, in light of other work I’ve done, doing Seurat-appropriate single-cell analysis outside of Seurat, here, we are now going to produce the same plot above using ggplot2 and another package called patchwork. We will stick with our so-called cytobank color scheme. The solution is below. Feel free to steal this code and use it however you’d like.

# Load necessary libraries
library(ggplot2)
library(patchwork)

# Define Cytobank-like color gradient
cytobank_colors <- colorRampPalette(c("blue", "cyan", "green", "yellow", "orange", "red"))

# Extract embedding and gene expression data
embedding <- Embeddings(pbmc, reduction = "umap")   
gene_data <- FetchData(pbmc, vars = genes)  
data <- cbind(as.data.frame(embedding), gene_data)

# Rename columns 
colnames(data) <- c("UMAP1", "UMAP2", genes)

# Initialize an empty list to store plots
plot_list <- list()

# Loop through each gene and create a plot
for (gene in genes) {
  p <- ggplot(data, aes(x = UMAP1, y = UMAP2, color = .data[[gene]])) +
    geom_point(size = 1, alpha = 0.8) +
    scale_color_gradientn(colors = cytobank_colors(100)) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 10),  # Smaller title font
      axis.title = element_text(size = 8),  # Smaller axis title font
      axis.text = element_text(size = 6),   # Smaller axis tick font
      legend.title = element_text(size = 8), # Smaller legend title font
      legend.text = element_text(size = 6)   # Smaller legend text font
    ) +
    guides(
      color = guide_colorbar(barwidth = 0.5, barheight = 4) # Adjust legend dimensions
    ) +
    labs(
      title = paste(gene),
      x = "UMAP1",
      y = "UMAP2",
      color = ""
    )
  plot_list[[gene]] <- p
}
# Combine all plots using patchwork
final_plot <- wrap_plots(plot_list, ncol = 3)  # Adjust ncol as needed
print(final_plot)