#Institute: Department of Immunology, Harvard Medical School
#Email: fah4zan@gmail.com, Fahd.Alhamdan@childrens.harvard.edu
#Abstact: Single-cell RNA sequencing (scRNA-Seq) has significantly enhanced our ability to study cellular heterogeneity and identify novel cell populations. However, a major challenge remains in bridging transcriptomic insights with protein-level validation using flow cytometry. To address this, I developed scFACS, an R package that transforms scRNA-Seq data into FACS-like plots, enabling researchers to visualize surface marker gene expression at the single-cell level before performing experimental validation. scFACS integrates seamlessly with existing single-cell analysis pipelines, including Seurat, to provide a user-friendly and interpretable approach for marker validation. The package allows users to gate cells based on marker expression, mimicking the workflow of traditional flow cytometry, and facilitates the identification of candidate markers for fluorescence-activated cell sorting (FACS). Additionally, scFACS enables the quantification of specific cell types, defined by known or newly identified surface markers, across different experimental conditions, providing insights into cellular dynamics and immune responses. In this example, I applied scFACS to single-cell RNA-Seq datasets derived from peripheral blood neutrophils collected from adults with sepsis — including both Survivors and Non-Survivors — as well as healthy controls (HC). For visualization, I used ITGAM (CD11b) as a canonical neutrophil marker and CD55 as a sepsis-associated marker. By bridging the gap between transcriptomics and proteomics, scFACS serves as a critical resource for researchers in immunology, cancer biology, and regenerative medicine, enabling efficient marker discovery and validation for flow cytometry-based studies.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 'SeuratObject' was built under R 4.4.0 but the current version is
## 4.4.1; it is recomended that you reinstall 'SeuratObject' as the ABI
## for R may have changed
## 'SeuratObject' was built with package 'Matrix' 1.7.0 but the current
## version is 1.7.3; it is recomended that you reinstall 'SeuratObject' as
## the ABI for 'Matrix' may have changed
##
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
##
## intersect, t
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(patchwork)
library(ggplot2)
library(RColorBrewer)
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
library(scFACS)
#load data
Adult.Neu <- readRDS("/Users/fahdalhamdan/Downloads/Adult.h.N.rds")
# Compute module scores
Adult.Neu <- AddModuleScore(Adult.Neu, features = list("ITGAM"), name = "x")
Adult.Neu <- AddModuleScore(Adult.Neu, features = list("CD55"), name = "y")
# Extract module scores for plotting
data <- FetchData(Adult.Neu, vars = c("x1", "y1", "condition"))
# Subset the two condition
data.HC <-data[data$condition == "HC", ]
data.S <-data[data$condition == "S", ]
data.NS <-data[data$condition == "NS", ]
# assign the quadrants
data.HC$quadrant <- with(data.HC, ifelse(x1 <= 0 & y1 > 0, "Q1",
ifelse(x1 > 0 & y1 > 0, "Q2",
ifelse(x1 < 0 & y1 < 0, "Q3", "Q4"))))
data.S$quadrant <- with(data.S, ifelse(x1 <= 0 & y1 > 0, "Q1",
ifelse(x1 > 0 & y1 > 0, "Q2",
ifelse(x1 < 0 & y1 < 0, "Q3", "Q4"))))
data.NS$quadrant <- with(data.NS, ifelse(x1 <= 0 & y1 > 0, "Q1",
ifelse(x1 > 0 & y1 > 0, "Q2",
ifelse(x1 < 0 & y1 < 0, "Q3", "Q4"))))
# Log transformation
data.HC$x1 <- log1p(data.HC$x1)
## Warning in log1p(data.HC$x1): NaNs produced
data.HC$y1 <- log1p(data.HC$y1)
## Warning in log1p(data.HC$y1): NaNs produced
data.S$x1 <- log1p(data.S$x1)
## Warning in log1p(data.S$x1): NaNs produced
data.S$y1 <- log1p(data.S$y1)
## Warning in log1p(data.S$y1): NaNs produced
data.NS$x1 <- log1p(data.NS$x1)
## Warning in log1p(data.NS$x1): NaNs produced
data.NS$y1 <- log1p(data.NS$y1)
## Warning in log1p(data.NS$y1): NaNs produced
# Calculate percentages in each quadrant
pct1 <- prop.table(table(data.HC$quadrant)) * 100
pct2 <- prop.table(table(data.S$quadrant)) * 100
pct3 <- prop.table(table(data.NS$quadrant)) * 100
print(pct1)
##
## Q1 Q2 Q3 Q4
## 21.590909 4.545455 64.772727 9.090909
print(pct2)
##
## Q1 Q2 Q3 Q4
## 42.937853 9.604520 38.418079 9.039548
print(pct3)
##
## Q1 Q2 Q3 Q4
## 40.577250 23.938879 25.806452 9.677419
# Create color palettes for each quadrant
colors <- colorRampPalette(rev(brewer.pal(11, 'Spectral')))
# Determine the limits for axes to perfectly fit the densities and positon the Qs
x_limits <- c(-4,3)
y_limits <- c(-4,4)
# Combine plots into a single density plot centered at 0 with percentage annotations around the center
p1 <- scFACS(data.HC) +
labs(x = "CD11b/ITGAM", y = "CD55", title = "HC", fill = "Density") +
# Adjust annotation positions near the center
annotate("text", x = 2, y = 4, label = paste(round(pct1["Q2"], 1), "%")) +
annotate("text", x = -2, y = 4, label = paste(round(pct1["Q1"], 1), "%")) +
annotate("text", x = -2, y = -3.9, label = paste(round(pct1["Q3"], 1), "%")) +
annotate("text", x = 2, y = -3.9, label = paste(round(pct1["Q4"], 1), "%")) +
xlim(x_limits) + # Adjust x-axis limits to fit the densities perfectly
ylim(y_limits) + # Adjust y-axis limits to fit the densities perfectly
theme_minimal() +
theme(axis.line = element_line(size = 1),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p2 <- scFACS(data.S) +
labs(x = "CD11b/ITGAM", y = "CD55", title = "Survivors", fill = "Density") +
# Adjust annotation positions near the center
annotate("text", x = 2, y = 4, label = paste(round(pct2["Q2"], 1), "%")) +
annotate("text", x = -2, y = 4, label = paste(round(pct2["Q1"], 1), "%")) +
annotate("text", x = -2, y = -3.9, label = paste(round(pct2["Q3"], 1), "%")) +
annotate("text", x = 2, y = -3.9, label = paste(round(pct2["Q4"], 1), "%")) +
xlim(x_limits) + # Adjust x-axis limits to fit the densities perfectly
ylim(y_limits) + # Adjust y-axis limits to fit the densities perfectly
theme_minimal() +
theme(axis.line = element_line(size = 1),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
p3 <- scFACS(data.NS) +
labs(x = "CD11b/ITGAM", y = "CD55", title = "Non-Survivors", fill = "Density") +
# Adjust annotation positions near the center
annotate("text", x = 2, y = 4, label = paste(round(pct3["Q2"], 1), "%")) +
annotate("text", x = -2, y = 4, label = paste(round(pct3["Q1"], 1), "%")) +
annotate("text", x = -2, y = -3.9, label = paste(round(pct3["Q3"], 1), "%")) +
annotate("text", x = 2, y = -3.9, label = paste(round(pct3["Q4"], 1), "%")) +
xlim(x_limits) + # Adjust x-axis limits to fit the densities perfectly
ylim(y_limits) + # Adjust y-axis limits to fit the densities perfectly
theme_minimal() +
theme(axis.line = element_line(size = 1),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank())
# Combine your ggplots
p1 + p2 + p3
