Unsupervised flow cytometry with R

#### Script used to generate figures and tables for the paper:
#"Viral and immune dynamics of HPV genital infections" currently in review (2024) in Plos Biology

# Original contributors:
# Nicolas Tessandier
# Baptiste Elie
# Samuel Alizon
# For any questions, please contact Nicolas Tessandier (email above or at https://tessandier.fr/contact/)
# Associated data can be found at: https://doi.org/10.57745/KJGOYZ
#####################
# NO COMMERCIAL USE #
# CC BY-NC-SA 4.0   #
#####################


##### Original script
suppressPackageStartupMessages(library('tidyverse'))
suppressPackageStartupMessages(library('magrittr'))
suppressPackageStartupMessages(library('flowCore'))
suppressPackageStartupMessages(library('FlowSOM'))
suppressPackageStartupMessages(library('CATALYST'))
suppressPackageStartupMessages(library('diffcyt'))
suppressPackageStartupMessages(library('cytofWorkflow'))
suppressPackageStartupMessages(library('ComplexHeatmap'))
suppressPackageStartupMessages(library('RColorBrewer'))
suppressPackageStartupMessages(library('ggpubr'))
suppressPackageStartupMessages(library('xtable'))
suppressPackageStartupMessages(library('viridis'))
suppressPackageStartupMessages(library('ggridges'))

# Load metadata
rm(list=ls())
experiment_info <- read_csv("./tables/experiment_info_PlosBio_R2.csv")
panel <- read_csv("./tables/panel_info_PlosBio_R2.csv")

# Load files
fcsfiles <- list.files(path = "./fcs", pattern = "\\.fcs$", full.names = TRUE)

# Load files in flowSet
d_flowSet <- read.flowSet(fcsfiles, transformation = FALSE, truncate_max_range = FALSE)

# Load fcs into sce

sce <- prepData(d_flowSet, panel, experiment_info,features = panel$fcs_colname, FACS = TRUE,
    md_cols = list(file = "file_name", id = "sample_id",
        factors = c("condition", "vaccine", "patient_id")))

# Set seed
set.seed(42)

# Cluster all cells
sce <- cluster(sce, features = rownames(sce)[c(1, 3, 6:7, 10, 13:15)],
xdim = 10, ydim = 10, maxK = 20, seed = 42, verbose = TRUE)

# Expression check (Figure just used to check CD45neg clusters,
# before restricting the dataset to CD45+ cells, not used in the paper)
pdf(paste0("./output/Figure0.pdf"), width = 16, height = 10)
u <- plotClusterHeatmap(sce,
hm2 = NULL, k = "meta20", m = NULL,
cluster_anno = TRUE, draw_freqs = TRUE)
plot(u)
dev.off()

# Remove CD45 negative cells

merging_manual_all_cells <- tibble(original_cluster = c(1:20),
new_cluster = c(
"Non-Leukocytes" ,# 1
"Non-Leukocytes" ,# 2
"Non-Leukocytes" ,# 3
"Non-Leukocytes" ,# 4
"Non-Leukocytes" ,# 5
"Non-Leukocytes" ,# 6
"Non-Leukocytes" ,# 7
"CD16+ leukocytes" ,# 8
"CD45low cells" ,# 9
"CD16+ granulocytes" ,# 10
"CD16+ granulocytes" ,# 11
"CD45low cells" ,# 12
"CD16+ granulocytes" ,# 13
"CD45+ CD3- CD16- cells (1)" ,# 14
"TCRgd cells" ,# 15
"CD16+ granulocytes" ,# 16
"CD16+ granulocytes" ,# 17
"CD45+ CD3- CD16- cells (2)" ,# 18
"CD4 T cells" ,# 19
"CD8 T cells" # 20
))

merging_manual_all_cells$new_cluster <- factor(merging_manual_all_cells$new_cluster,
levels = c(
"Non-Leukocytes" ,# 1
"CD16+ leukocytes" ,# 8
"CD45low cells" ,# 9
"CD16+ granulocytes" ,# 10
"CD45+ CD3- CD16- cells (1)" ,# 14
"TCRgd cells" ,# 15
"CD45+ CD3- CD16- cells (2)" ,# 18
"CD4 T cells" ,# 19
"CD8 T cells" # 20
))
merging_manual_all_cells

### apply manual merging ###
sce <- mergeClusters(sce, k = "meta20", table = merging_manual_all_cells, id = "merging1", overwrite=TRUE)

sub <- filterSCE(sce, !cluster_id %in% c(1, 2, 3, 4, 5, 6, 7) , k = "meta20")

# Re-clustering the CD45-restricted subset
set.seed(42)
sub <- cluster(sub, features = rownames(sce)[c(1, 3, 6:7, 10, 13:15)],
xdim = 10, ydim = 10, maxK = 20, seed = 42)

#
# Annotating clusters check (Figure just used to check manual annotation of clusters, not used in the paper)
pdf(paste0("./output/Figure02.pdf"), width = 16, height = 10)
u <- plotClusterHeatmap(sub,
hm2 = NULL, k = "meta20", m = NULL,
cluster_anno = TRUE, draw_freqs = TRUE)
plot(u)
dev.off()

merging_manual_leuko_only <- tibble(original_cluster = c(1:20),
new_cluster = c(
"I (CD16+ granulocytes)" , # 1
"I (CD16+ granulocytes)" , # 2
"III (CD16+ leukocytes)" , # 3
"II (CD16+ granulocytes)" , # 4
"I (CD16+ granulocytes)" , # 5
"XI (CD8 T cells)" , # 6
"VII (CD45low cells)" , # 7
"II (CD16+ granulocytes)" , # 8
"VIII (CD45low cells)" , # 9
"IX (TCRgd cells)" , # 10
"IV (CD16+ leukocytes)" , # 11
"I (CD16+ granulocytes)" , # 12
"IX (TCRgd cells)" , # 13
"V (CD45+ CD16- CD3- cells)" , # 14
"VIII (CD45low cells)" , # 15
"VIII (CD45low cells)" , # 16
"X (CD4 T cells)" , # 17
"V (CD45+ CD16- CD3- cells)" , # 18
"VI (CD45+ CD16- CD3- cells)" , # 19
"II (CD16+ granulocytes)" # 20
))

merging_manual_leuko_only$new_cluster <- factor(merging_manual_leuko_only$new_cluster,
levels = c(
"I (CD16+ granulocytes)" , # 1
"II (CD16+ granulocytes)" , # 2
"III (CD16+ leukocytes)" , # 3
"IV (CD16+ leukocytes)" , # 8
"V (CD45+ CD16- CD3- cells)" , # 14
"VI (CD45+ CD16- CD3- cells)" , # 7
"VII (CD45low cells)" , # 9
"VIII (CD45low cells)" , # 15
"IX (TCRgd cells)" , # 13
"X (CD4 T cells)" , # 17
"XI (CD8 T cells)"  # 6
))
merging_manual_leuko_only

### Merging for dataset restricted to leukocytes
sub <- mergeClusters(sub, k = "meta20", table = merging_manual_leuko_only, id = "merging1", overwrite=TRUE)

# Color palette for all the FCM figures
q11 <- CATALYST:::.cluster_cols[c( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)]

### UMAP
set.seed(42)
sub <- runDR(sub, "UMAP", cells = 150, features = rownames(sub)[c(1,3,7,6,14,13,15,10,8,11)])

## UMAP by condition
x <- plotDR(sub, "UMAP", color_by = "merging1", k_pal =  q11) +
facet_wrap("condition", labeller = labeller(vs = label_both, am = label_value) ) +
theme(legend.text=element_text(size=13),legend.position = "bottom") +
guides(colour = guide_legend(nrow =3,override.aes = list(size=3))) +
theme(strip.text.x = element_text(size = 20, face = "bold"),
strip.text.y = element_text(size = 20,face = "bold")) +
theme(legend.title=element_blank())

### Differential Abundance analysis ###
fdr_cutoff <- 0.1

ei <- metadata(sub)$experiment_info
da_formula <- createFormula(ei, cols_fixed = c("condition", "vaccine"), cols_random = "patient_id")
ds_design <- createDesignMatrix(ei, cols_design = c("condition", "vaccine"))
contrast.focal <- createContrast(c(0, 1, 0))
contrast.vaccine <- createContrast(c(0, 0, 1))

# Define the variables
da_res1_leuko_only.edgeR.focal <- diffcyt(
       sub,
       experiment_info = experiment_info,
       marker_info = panel,
       design = ds_design,
       formula = da_formula,
       contrast = contrast.focal,
       analysis_type = c("DA"),
       method_DA = c("diffcyt-DA-edgeR"),
       clustering_to_use = "merging1",
       plot = TRUE,
       path = "./output",
       verbose = TRUE
     )

da_res1_leuko_only.edgeR.vaccine <- diffcyt(
       sub,
       experiment_info = experiment_info,
       marker_info = panel,
       design = ds_design,
       formula = da_formula,
       contrast = contrast.vaccine,
       analysis_type = c("DA"),
       method_DA = c("diffcyt-DA-edgeR"),
       clustering_to_use = "merging1",
       plot = TRUE,
       path = "./output",
       verbose = TRUE
     )

plasma_pal <- c(viridis::viridis(n = 10, alpha = 0.6))
da_table1 <- as_tibble(topTable(da_res1_leuko_only.edgeR.focal, all = TRUE,
show_props = TRUE, format_vals = TRUE, digits = 3))

# Figure 2b
rowData(da_res1_leuko_only.edgeR.focal$res) %>% as_tibble() %>%
dplyr::select(cluster_id, logFC, p_val, p_adj) %>% mutate(FoldChange = 2^(logFC)) %>%
dplyr::select(-c(logFC)) %>% arrange(FoldChange, p_adj) %>% relocate(FoldChange, .before = p_val) -> res.edgeR

da_table1 %>%
	dplyr::select(-c(p_val)) %>%
	left_join(res.edgeR, by = c("cluster_id")) %>%
	dplyr::select(!c(p_adj.x, p_val)) %>%
	pivot_longer(!c(cluster_id, p_adj.y, FoldChange), names_to = "sample", values_to = "value") %>%
	mutate(sample_id = str_sub(sample, 7, 14)) %>%
	group_by(sample_id) %>%
	left_join(experiment_info) %>%
 ggplot(aes(fill= fct_rev(condition), y=value, x=cluster_id, label = FoldChange)) +
    geom_boxplot(outlier.shape = NA, reverse=TRUE) +
    scale_fill_manual( values = plasma_pal[c(10,1)], name = "") +
    scale_color_manual(values = plasma_pal[c(10,1)], name = "") +
    geom_point(colour = "black",position=position_jitterdodge(jitter.width = 0.15), alpha=0.5, size = 0.5) +
    #geom_text(aes(label = paste0("FC: ",round(FoldChange, 2),
    "\np.adj = ", round(p_adj.y, 4)), y=97, x = cluster_id) ) +
coord_flip() +
guides(fill = guide_legend(reverse = TRUE)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 0, hjust = 1, vjust = 1)) +
theme(legend.text=element_text(size=13, face = "bold"),legend.position = "bottom") +
theme(axis.text.x = element_text(face = "bold", size = 13)) +
theme(axis.text.y = element_text(face = "bold", size = 13))  +
scale_x_discrete(labels= c("1","2","3","4","5","6","7","8","9","10","11")) +
    xlab("") +
    ylim(c(0, 105)) +
    ylab("Frequency (%)") -> fig.2b

# Table for figure 2C
rowData(da_res1_leuko_only.edgeR.focal$res) %>% as_tibble() %>%
dplyr::select(cluster_id, logFC, p_val, p_adj) %>% mutate(FoldChange = 2^(logFC))
%>% dplyr::select(-c(logFC)) %>% arrange(FoldChange, p_adj)
%>% relocate(FoldChange, .before = p_val) %>% mutate(across(where(is.numeric), round, 3))


# Empty plot for figure 2C, the grid is generated manually through inkscape
# (and a lot of the figure is adapted that way, the legend of the UMAP for figure 2A, the Y axis for figure 2B)
bb <- ggplot(mtcars, aes(x = wt, y = mpg)) + geom_blank() + theme_void()
line2 <- ggarrange(x, widths = c(1),nrow = 1, labels =  c("A"),
font.label = list(size = 18, color = "black", face = "bold", family = NULL))
line3 <- ggarrange(fig.2b,bb, widths = c(2, 1),nrow = 1,
labels =  c("B", "C"),font.label = list(size = 18, color = "black", face = "bold", family = NULL))
figure2 <- ggarrange(line2, line3, ncol = 1)
ggsave(paste0("output/figure1_PlosBio_Tessandier_Elie.pdf"), plot = figure2,
       width = 15,
       height = 15,
       limitsize = FALSE,
       units = "in",
       dpi = 300)

# save FCM environement generated with figure1_PlosBio_Tessandier_Elie.r
save(file = "./RData/FCM_analysis_Plos_Bio_Tessandier_Elie.RData")