In jonesnoaht/BTIPFlow: What the Package Does (This Package Assists with Flow Analysis)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(tidyverse)
library(BTIPFlow)
library(flowCore)
library(flowViz)
library(flowWorkspace)
library(ggcyto)
library(cytotidyr)

get_fcs_resultsQC(QC_folder = "./RG Microbiome -APCs panel/resultsQC/", 
                  raw_folder = "../20200810NJ-ACEFlow/good/") -> cs

data(cs)
cs[[1]]@exprs[cs[[1]] %in% good_events]
blue[grepl("Comp", sampleNames(blue))]

data(cs)
gs <- GatingSet(cs)
sampleNames(gs)[6] <- "APCs_CD11C FMO_QC.fcs"
sampleNames(gs)[9] <- "APCs_CD45 FMO_QC.fcs"
markernames(gs[[1:24]])[[4]] <- "CD80"
markernames(gs[[1:24]])[[7]] <- "Live dead"

pData(gs) <- cbind(pData(gs),row.names(pData(gs)))
colnames(pData(gs))[2] <- "tube names"
for (i in 1:length(gs)) {
  pData(gs[[i]])$`tube names` <- keyword(gs[[i]])$`TUBE NAME`
}

good_events <- rectangleGate(filterId = "Good Events", "remove_from_all" = c(0,9999))
bad_events <- rectangleGate(filterId = "Bad Events", "remove_from_all" = c(10000,19999))

gs_pop_add(gs, good_events)
gs_pop_add(gs, bad_events)

set_scatter_gate(gating_set = gs,
                 gate_coordinates = c(2.0e+04, 2.0e+04,
                                      1.4e+04, 1.5e+04,
                                      1.4e+04, 0.2e+04,
                                      3.0e+04, 0.2e+04,
                                      3.0e+04, 2.0e+04),
                 gate_name = "scatter")
set_scatter_gate(gating_set = gs, 
                 gate_coordinates = c(0.01e+05, 1.5e+05,
                                      0.01e+05, 1.0e+05,
                                      0.12e+05, 1.0e+05,
                                      0.12e+05, 1.5e+05),
                 gate_name = "beads"
                 )

set_scatter_gate(gating_set = gs, 
                 gate_coordinates = c(0.7e+05, 0.4e+05,
                                      0.7e+05, 0.01e+05,
                                      0.95e+05, 0.01e+05,
                                      0.95e+05, 0.4e+05),
                 gate_name = "FSC-singlet", 
                 parent = "scatter",
                 dimensions = list("FSC-W", "FSC-H")
                 )

set_scatter_gate(gating_set = gs, 
                 gate_coordinates = c(0.4e+05, 0.3e+05,
                                      0.4e+05, 0.01e+05,
                                      0.70e+05, 0.01e+05,
                                      0.70e+05, 0.3e+05),
                 gate_name = "SSC-singlet", 
                 parent = "FSC-singlet",
                 dimensions = list("SSC-W", "SSC-H")
                 )

make_gates_from_fmos(gating_set = gs, parent = "SSC-singlet", trim = 1)

pdf("gating_tree.pdf")
plot(gs)
dev.off()

nodes <- gs_get_pop_paths(gs, path = "auto")[2:4]
length(gs) %>% sqrt() %>% round() -> hw
hw <- hw*2
pdf("Root Scatter Plots.pdf", hw, hw)
ggcyto(gs, 
       aes(`FSC-A`, `SSC-A`), 
       subset = "root", 
       filter = marginalFilter) + 
  geom_bin2d(bins = 256) + 
  ggcyto_par_set(limits = "instrument") +
  geom_gate("scatter") + 
  geom_stats("scatter", type = "percent", adjust = 0.8) +
  labs(title = "Root Scatter Plots") +
  facet_wrap(~`tube names`)
dev.off()
pdf("FSC-singlet Plots.pdf", hw, hw)
ggcyto(gs, 
       aes(`FSC-W`, `FSC-H`), 
       subset = "scatter", 
       filter = marginalFilter) + 
  geom_bin2d(bins = 256) + 
  ggcyto_par_set(limits = "instrument") +
  geom_gate("FSC-singlet") + 
  geom_stats("FSC-singlet", type = "percent", adjust = 0.1) +
  labs(title = "FSC-singlet Plots") +
  facet_wrap(~`tube names`)
dev.off()
pdf("SSC-singlet Plots.pdf", hw, hw)
ggcyto(gs, 
       aes(`SSC-W`, `SSC-H`), 
       subset = "FSC-singlet", 
       filter = marginalFilter) + 
  geom_bin2d(bins = 256) + 
  ggcyto_par_set(limits = "instrument") +
  geom_gate("SSC-singlet") + 
  geom_stats("SSC-singlet", type = "percent", adjust = 0.1) +
  labs(title = "SSC-singlet Plots") +
  facet_wrap(~`tube names`)
dev.off()

pops <- gs_get_pop_paths(gs, path = "auto")
for (i in 7:13) {
  file_name <- paste(pops[i], ".pdf", sep = "")
  param <- (gs_pop_get_gate(gs[[1]], pops[i])[[1]] %>% parameters())[[1]]
  ggcyto(gs, 
         aes_(param, "SSC-A"), 
         subset = pops[6]) + 
    geom_bin2d(bins = 256) + 
    scale_x_flowjo_biexp(neg = 1, widthBasis = -20) +
    geom_gate(pops[i]) + 
    geom_stats(pops[i], type = "percent", adjust = 0.1) +
    labs(title = paste(pops[i], "Plots"), 
         subtitle = "Gated per FMO") +
    facet_wrap(~`tube names`)
  ggsave(file_name, width = hw, height = hw)
}

Now you can plot your figures

First we extract our data from our gatingset.

stats <- gs_pop_get_count_with_meta(gs)
stats$`group` <- ""
for (i in 1:length(stats$`tube names`)) {
  if (grepl("ABX", stats[i]$`tube names`)) {
    stats[i]$group <- "Antibiotic"
  } else {
    if (grepl("Ctrl", stats[i]$`tube names`) & !grepl("spleen", stats[i]$`tube names`)) {
      stats[i]$group <- "Control"
    }
  }
}

Next we refine our data and generate a new column with absolute bead counts.

abs <- data.frame()
stats %>% dplyr::filter(group != "") -> stats_good
stats_good$`tube name` %>% factor() %>% levels() -> levels_i
(stats_good %>% dplyr::filter(grepl("pos", Population)))$Population %>% factor() %>% levels() -> levels_j
for (i in levels_i) {
  for (j in levels_j) {
    x <- c(i, j, ((stats_good %>% dplyr::filter(`tube names` == i, `Population` == j) %>% select(`Count`))[[1]] / (stats_good %>% dplyr::filter(`tube names` == i, `Population` == "/beads") %>% select(`Count`))[[1]] * 1e+04))
    abs <- rbind(abs, x)
  }
}
abs <- cbind(abs, rep("", nrow(abs)))
colnames(abs) <- c("tube", "pop", "abs_count", "group")
abs$abs_count <- as.numeric(abs$abs_count)
for (i in 1:length(abs$`tube`)) {
  if (grepl("ABX", abs$`tube`[i])) {
    abs$group[i] <- "Antibiotic"
  } else {
    if (grepl("Ctrl", abs$`tube`[i])) {
      abs$group[i] <- "Control"
    }
  }
}

The default plot titles will be terrible unless we do something about the population names.

abs$pop <- sub("/scatter/FSC-singlet/SSC-singlet/Live dead neg/", "", abs$pop)
abs$pop <- sub("/scatter/FSC-singlet/SSC-singlet/", "", abs$pop)

Then we generate the plots.

ggplot(abs, aes(x = group, y = abs_count)) + 
  geom_jitter() + 
  ylim(0, 1.4*max(abs$abs_count)) +
  stat_summary(geom = "errorbar", 
                 width = .2, 
                 position = position_dodge(), 
                 color = 2, 
                 fun.data = mean_se
                 ) +
    ggsignif::geom_signif(comparisons = list(c("Control", "Antibiotic")),
      test = "t.test",
      y_position = max(abs$abs_count), 
      step_increase = .1
    ) + 
    ggpubr::theme_pubr() +
  facet_wrap(~`pop`)
ggsave("ctrl_vs_abx.pdf", width = 7)

jonesnoaht/BTIPFlow documentation built on Feb. 20, 2024, 8:45 a.m.