Dye Order & Distinguishing Power

Background

This is the flow cytometry data for miglog phase Cg stained with PI/PI+SYTO9/ SYTO+PI in different staining orders, and run through flow cytometry (details in ELN). The goal is to establish FungaLight’s distinguishing power with different staining protocols.

require(tidyverse)
require(flowCore)
require(flowClust)
require(openCyto)
require(ggcyto)
require(cowplot)
require(ggridges)

old <- theme_set(theme_minimal())

Import data

FCS files are stored in RDSS/user/flow cytometry, FCS is read and write into an input tsv table. The tsv file is avaliable from the Input folder.

# use relative path to make it easier for collaboration
data.path = "/Volumes/rdss_bhe2/User/Hanxi Tang/Flow Cytometry/20240109_0116 Dye order/"
dat0 <- read.flowSet(path = data.path, transformation = FALSE,  # the original values are already linearized. 
                     emptyValue = FALSE,  alter.names = TRUE,   # change parameter names to R format
                     column.pattern = ".H|FSC|SSC") # only load the height variables for the fluorescent parameters

Simplify the sample names

#require(PTXQC)
#source("../script/20220326-simplify-names-subroutine.R")
oriNames <- sampleNames(dat0)
tmp <- str_split(oriNames, pattern = "[ _]+", simplify = TRUE)[,c(1, 6, 7)]
colnames(tmp) <- c("Date", "Treatment", "Dye")
sample <- data.frame(tmp) %>% 
  mutate(
    Dye = case_match(
      Dye,
      "p.fcs" ~ "PI alone",
      "pf.fcs" ~ "PI first",
      "sf.fcs" ~ "SYTO9+PI",
      .default = NA
    )
  )
rownames(sample) <- oriNames
pData(dat0) <- sample
dat <- dat0[!is.na(sample$Dye)]
print(pData(dat) %>% as_tibble())

write_tsv(pData(dat), file = "../data/10.01-sample-list.tsv")

Basic gating

Outlier

The following gaphing steps are used to gate singlets by FSC and SSC values. Only singlets are included in analysis.

outlier.gate <- rectangleGate(filterId = "-outlier", "FSC.H" = c(1.2e5, 1e6), "SSC.H" = c(1e2, 1e6))
ggcyto(dat[c(2,20)], aes(x = FSC.H, y = SSC.H), subset = "root") +
  geom_hex(bins = 64) + geom_gate(outlier.gate) + facet_wrap(~name, ncol = 2) + ggcyto_par_set(limits = "instrument")

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Add gate to GS

# shorten the names of the samples for plotting
shortNames <- with(pData(dat), 
                   paste(substr(Date, 1, 4), paste0(Treatment, " mM")
                         , Dye, sep = "_"))
names(shortNames) <- sampleNames(dat)
gs <- GatingSet(dat) # create a GatingSet
# rename the GatingSet
if(all(sampleNames(gs) == names(shortNames)))
  pData(gs)$name <- shortNames
# repair the Date column in pData
pData(gs)$Date = paste0("0", pData(gs)$Date)
# add the outlier gate
gs_pop_add(gs, outlier.gate, parent = "root")

[1] 2

# compute
recompute(gs)

done!

Gate for singlets

scPars <- ggcyto_par_set(limits = list(x = c(0,1e6), y = c(30,300)))
ex <- Subset(dat[[1]], outlier.gate)
polygon <- matrix(c(1e5, 1e5, 1e6, 1e6, 
                    60, 105, 135,60), ncol = 2)
colnames(polygon) <- c("FSC.H", "FSC.W")
singlet.gate <- polygonGate(filterId = "singlet", .gate = polygon)
ggcyto(ex, aes(x = FSC.H, y = FSC.W)) + geom_hex(bins = 128) + geom_gate(singlet.gate) + geom_stats() + scPars

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Add this gate to the gatingSet

gs_pop_add(gs, singlet.gate, parent = "-outlier", name = "singlet")

[1] 3

recompute(gs)

done!

scPars <- ggcyto_par_set(limits = list(x = c(0,1e6), y = c(30,200)))
test <- sample(1:length(gs), 8)
ggcyto(gs[test], aes(x = FSC.H, y = FSC.W), subset = "-outlier") +
  geom_hex(bins = 128) + geom_gate("singlet") + facet_wrap(~name, ncol = 4) + 
  scPars# + theme(strip.text = element_blank())

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

FSC and SSC vs treatment and dye

FSC

mult_format <- function() {
     function(x) format(x/10000,digits = 2) 
}
for(date in c("011124", "011524", "011624")){
  p <- ggcyto(gs[pData(gs)$Date == date], aes(x = FSC.H), subset = "singlet") + 
    geom_density_ridges(aes(y = Treatment), fill = "forestgreen", alpha = 0.8) + 
    scale_x_continuous(labels = mult_format(), name = "FSC.H x 10000") +
    facet_wrap(~Dye) + labs(title = paste0("Date: ", date))
  print(p)
}

Mock treated sample stands out - the three H2O2 treated samples show a trend of decreasing FSC with increasing ROS Mock has a FSC in between 10 and 100 mM

FSC vs BL1.H (green staining)

for(date in c("011124", "011524", "011624")){
  p <- ggcyto(gs[pData(gs)$Date == date], aes(x = BL1.H, y = FSC.H), subset = "singlet") + 
    geom_hex(bins = 128) + 
    #geom_gate(live.gate) + 
    #geom_stats(location = "data", adjust = c(0.1, 0.8), digits = 2) + 
    labs(title = paste0("Date: ", date)) +
    facet_grid(Dye ~ Treatment) +
    scale_x_logicle() + #scale_y_logicle() + 
    theme(strip.text = element_text(size = rel(1.1), color = "steelblue", face = 2))
  print(p)
}

SSC

for(date in c("011124", "011524", "011624")){
  p <- ggcyto(gs[pData(gs)$Date == date], aes(x = SSC.H), subset = "singlet") + 
    geom_density_ridges(aes(y = Treatment), fill = "forestgreen", alpha = 0.8) + 
    scale_x_continuous(labels = mult_format(), name = "FSC.H x 10000") +
    facet_wrap(~Dye) + labs(title = paste0("Date: ", date))
  print(p)
}

SSC (complexity of cell, roughly) decreases with increasing H2O2 concentration Mock again show a different pattern.

Gate for scoring

Live = red negative

This is the live population with red fluorescence below 10^2.2

Important: I changed the Red top boundary to 10^2.3 from 10^2.2

Main takeaway

1. From 10mM to 1M treated samples, %live decreases and the change is much more obvious in the dual stain (PI first or simultaneous) than PI alone

2. The mock treated sample has a visible high red high green (not as high and prominent as the 100 mM treated) than 10 mM treated. It doesn’t follow the trend above.

#scPars <- ggcyto_par_set(limits = list(x = c(0,1e6), y = c(0,1e5)))
polygon <- matrix(c(0, 10^5, 10^5, 0,
                    10^2.3, 10^2.3, 0, 0), ncol = 2)
colnames(polygon) <- c("BL1.H", "BL3.H")
live.gate <- polygonGate(filterId = "live", .gate = polygon)
for(date in c("011124", "011524", "011624")){
  p <- ggcyto(gs[pData(gs)$Date == date], aes(x = BL1.H, y = BL3.H), subset = "singlet") + 
    geom_hex(bins = 128) + 
    geom_gate(live.gate) + 
    geom_stats(location = "data", adjust = c(0.1, 0.8), digits = 2) + 
    labs(title = paste0("Date: ", date)) +
    facet_grid(Dye ~ Treatment) +
    scale_x_logicle() + scale_y_logicle() + 
    theme(strip.text = element_text(size = rel(1.1), color = "steelblue", face = 2))
  print(p)
}

Notice the higher percentage of events in the live gate in 10 mM vs mock.

Add this gate to the gatingSet

Intermediate = high G/R

This is the intermediate population with high Red and Green. We hypothesize that they represent the

Important: I changed the bottom boundary from 10^2.2 to 10^2.3, and the left boundary from 10^3 to 2*10^3

Main takeaway

1. Different from the %live gate, the intermediate population is most obvious in the 100 mM treated sample, and absent in both the 10 mM and 1M treated ones. This suggests that this population could be useful for distinguishing stress levels below the lethal concentration but above the tolerance levels of the species - in this case, < 1M and > 10 mM.

2. The mock treated sample showed an intermediate population. Was unexpected.

#scPars <- ggcyto_par_set(limits = list(x = c(0,1e6), y = c(0,1e6)))
#ex <- Subset(dat.f1g[[9]], singlet.gate)
polygon <- matrix(c(10^5, 2*10^3, 10^5,
                    10^2.3, 10^2.3, 10^4), ncol = 2)
colnames(polygon) <- c("BL1.H", "BL3.H")
inter.gate <- polygonGate(filterId = "inter", .gate = polygon)
for(date in c("011124", "011524", "011624")){
  p <- ggcyto(gs[pData(gs)$Date == date], aes(x = BL1.H, y = BL3.H), subset = "singlet") + 
    geom_hex(bins = 128) + 
    geom_gate(inter.gate) + 
    geom_stats(location = "data",  adjust = c(0.1, 0.8), digits = 2) + 
    labs(title = paste0("Date: ", date)) +
    facet_grid(Dye ~ Treatment) +
    scale_x_logicle() + scale_y_logicle() + 
    theme(strip.text = element_text(size = rel(1.1), color = "steelblue", face = 2))
  print(p)
}

Add this gate to the gatingSet

Dead = High red low green

We don’t need to set the dead gate because all we want is the percent of events in that gate, not the MFI for either channel. The percentage can be calculated by subtracting the percent of events in the previous two gates from 100%.

Exported gate stats

% of events in each of the three gates

gated_stats <- gs_pop_get_stats(gs, type = "count") %>% 
  as_tibble() %>% 
  mutate(pop = gsub(".*/", "", pop), pop = gsub("-outlier", "cells", pop)) %>% 
  pivot_wider(names_from = pop, names_prefix = "n_", values_from = count) %>% 
  mutate(
    #p_singlet = n_singlet / n_cells,
    p_live = n_live / n_singlet,
    p_intermediate = n_intermediate / n_singlet,
    p_dead = 1 - p_live - p_intermediate
  )
#write_tsv(gated_data, file = "../Input/20240124-fig-2-Dye-Order-gated-ht.tsv")
#print(gated_data)

MFI in the live and intermediate gates

popMFI <- gs_pop_get_stats(gs, nodes = c("live", "intermediate"), type = pop.MFI) %>% 
  as_tibble() %>% 
  select(sample, pop, FSC = `FSC-H`, BL1 = `Fungalight Green-H`) %>% 
  pivot_wider(names_from = pop, values_from = c(FSC, BL1))

Combine the data and meta data

final <- select(pData(dat), sample = name, everything()) %>% 
  left_join(select(gated_stats, -starts_with("n_")),  by = "sample") %>% 
  left_join(popMFI, by = "sample")

write_tsv(final, file = "../data/12.01-gated-output-20240117.tsv")
print(final %>% select(-sample))

Exploratory comparison between PI alone and combo

PI alone

p1 <- final %>% 
  dplyr::filter(Dye == "PI alone") %>% 
  ggplot(aes(x = Treatment, y = 1-p_live)) + 
  geom_point(size = 2, shape = 21, position = position_jitter(0.05)) + 
  #stat_summary(fun = "mean", geom = "point", color = "red") +
  labs(x = bquote(H[2]*O[2]~(mM)), y = "1 - %Live") +
  theme_cowplot()

inset <- final %>% 
  dplyr::filter(Dye == "PI alone", Treatment != "1000") %>% 
  ggplot(aes(x = Treatment, y = 1-p_live)) + 
  geom_bar(stat = "summary", fun = "mean", fill = "skyblue2", alpha = 0.7) +
  geom_point(size = 2, position = position_jitter(0.05)) + 
  #stat_summary(fun = "mean", geom = "point", size = 15, shape = "-", color = "red") +
  labs(x = bquote(H[2]*O[2]~(mM)), y = "1 - %Live") +
  theme_minimal_hgrid(11)

ggdraw(p1 + labs(subtitle = "PI alone, doesn't disguish 10 vs 100 mM")) + 
  draw_plot(inset, .2, .8, .4, .5, vjust = 1) +
  draw_plot_label(c("A", "B"), c(0, 0.18), c(1, 0.85))

ggsave("../output/12.01-PI-alone-one-minus-perc-live.png", width = 5, height = 4)

PI+SYTO9

p2 <- final %>% 
  dplyr::filter(Dye != "PI alone", Treatment != "0") %>% 
  ggplot(aes(x = Treatment, y = 1-p_live)) + 
  geom_bar(stat = "summary", fun = "mean", fill = "gray", alpha = 0.7) +
  geom_point(size = 2, position = position_jitter(0.05)) + 
  #stat_summary(fun = "mean", geom = "point", shape = "-", 
  #             size = 15, color = "red", position = position_nudge(x = 0.2)) +
  labs(x = bquote(H[2]*O[2]~(mM)), y = "1 - %Live") +
  facet_wrap(~Dye) +
  theme_cowplot() +
  theme(strip.text = element_text(size = rel(1.2), face = 2),
        strip.background = element_blank())

p2
ggsave("../output/12.01-PI+SYTO9-one-minus-plive.png")

Saving 7.29 x 4.51 in image

Plot % Dead and Compare

```r
pf1g <- gated_data %>%
  dplyr::filter(!(Dye %in% c(\PI\))) %>% 
  dplyr::filter(!(Treatment %in% c(\1000 mM\)))  %>%
  #mutate(Dilution = factor(Dilution, levels = c(\250 x\, \500 x\, \1000 x\, \2000 x\))) %>%
  ggplot(aes(x = Treatment, y = percent_Live)) +
  facet_wrap(~ Dye, scale = \free_y\) +
  geom_boxplot() +
  geom_point(stroke = 1, size = 2, position = position_jitter(0.2)) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()
pf1g +
  labs(x = \Hydrogen Peroxide Treatment Concentration\, y = \Gated Live Cells\) +
  panel_border(color = \black\, size = 1.5) +
  theme(axis.line = element_blank(),
        strip.background = element_blank(),
        axis.text.x = element_text(size = 12, face = \bold\, color = \black\),  # Adjust size and face (bold)
        axis.text.y = element_text(size = 16, face = \bold\, color = \black\),
        axis.title.x = element_text(size = 16, face = \bold\),
        axis.title.y = element_text(size = 16, face = \bold\)
        )+ 
  theme(strip.text = element_text(size = 16, face = \bold\)) 
#ggsave(\../output/20240104-PI-Dilutions-Gated_RP-Plot.png\, width = 7, height = 5)

```