Flow Cytometry Analysis: Complete Guide & Free R Tutorial

# Essential packages for flow cytometry analysis
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Core flow cytometry packages
BiocManager::install(c(
    "flowCore",           # Core flow cytometry data structures
    "flowWorkspace",      # Workspace management
    "ggcyto",            # Grammar of graphics for cytometry
    "openCyto",          # Automated gating framework
    "CytoML"             # Import/export gating sets
))

# Statistical and visualization packages
install.packages(c(
    "ggplot2",           # Advanced plotting
    "dplyr",             # Data manipulation
    "reshape2",          # Data reshaping
    "RColorBrewer",      # Color palettes
    "gridExtra",         # Multi-panel plots
    "scales"             # Plot scaling functions
))

# Load essential libraries
library(flowCore)
library(ggcyto)
library(dplyr)
library(ggplot2)

cat("✅ Flow cytometry analysis environment ready!\n")

# Import single FCS file for flow cytometry analysis
flow_data <- read.FCS("sample_data.fcs", 
                     transformation = FALSE,
                     truncate_max_range = FALSE)

# Examine data structure and parameters
cat("📊 Flow Cytometry Data Summary:\n")
cat("Cells analyzed:", nrow(flow_data), "\n")
cat("Parameters measured:", ncol(flow_data), "\n")

# Display parameter information
parameter_info <- flow_data@parameters@data
print(parameter_info[, c("name", "desc", "range", "minRange", "maxRange")])

# Quality control checks
exprs_data <- exprs(flow_data)

# Check for anomalies
negative_values <- sum(exprs_data < 0)
zero_values <- sum(exprs_data == 0)
outliers <- sum(exprs_data > quantile(exprs_data, 0.999, na.rm = TRUE))

cat("\n🔍 Quality Control Results:\n")
cat("Negative values detected:", negative_values, "\n")
cat("Zero values:", zero_values, "\n")
cat("Potential outliers:", outliers, "\n")

# Time-based quality assessment
if("Time" %in% colnames(exprs_data)) {
    time_range <- range(exprs_data[,"Time"])
    acquisition_time <- (time_range[2] - time_range[1]) / 100 # centiseconds to seconds
    cat("Acquisition time:", round(acquisition_time, 2), "seconds\n")
}

# Create transformation functions for flow cytometry analysis
# Logarithmic transformation for fluorescence channels
log_transform <- logTransform(transformationId = "log10",
                             logbase = 10, r = 1, d = 1)

# Asinh transformation (alternative for low signal data)
asinh_transform <- arcsinhTransform(transformationId = "asinh",
                                   a = 1, b = 1, c = 1)

# Identify fluorescence channels (exclude FSC, SSC, Time)
fluor_channels <- colnames(flow_data)[!grepl("FSC|SSC|Time", 
                                           colnames(flow_data))]

# Create transformation list
transform_list <- transformList(fluor_channels, log_transform)

# Apply transformation
flow_transformed <- transform(flow_data, transform_list)

# Compare before and after transformation
par(mfrow = c(1, 2))
hist(exprs(flow_data)[, fluor_channels[1]], 
     main = "Before Transformation", 
     xlab = fluor_channels[1],
     breaks = 50, col = "lightblue")

hist(exprs(flow_transformed)[, fluor_channels[1]], 
     main = "After Log Transformation", 
     xlab = paste("Log10", fluor_channels[1]),
     breaks = 50, col = "lightgreen")

cat("✅ Data transformation completed successfully!\n")

# Compensation for flow cytometry analysis
library(flowStats)

# Load compensation controls (single-stained samples)
comp_files <- list.files(pattern = "comp.*\\.fcs$", full.names = TRUE)

# Create compensation matrix from single-stained controls
if(length(comp_files) > 0) {
    # Load compensation controls
    comp_samples <- lapply(comp_files, read.FCS)
    names(comp_samples) <- gsub("\\.fcs$", "", basename(comp_files))
    
    # Calculate compensation matrix
    comp_matrix <- spillover(comp_samples[[1]], 
                           unstained = comp_samples[["unstained"]],
                           patt = "FITC|PE|APC|PerCP")
    
    # Display compensation matrix
    print("🔧 Calculated Compensation Matrix:")
    print(comp_matrix)
    
    # Apply compensation
    flow_compensated <- compensate(flow_transformed, comp_matrix)
    
} else {
    # Use pre-calculated compensation if available
    comp_matrix <- flow_data@description$SPILL
    
    if(!is.null(comp_matrix)) {
        flow_compensated <- compensate(flow_transformed, comp_matrix)
        cat("✅ Applied pre-calculated compensation matrix\n")
    } else {
        flow_compensated <- flow_transformed
        cat("⚠️ No compensation matrix found - proceeding without compensation\n")
    }
}

# Visualize compensation effect
library(ggcyto)
if(length(fluor_channels) >= 2) {
    p1 <- ggcyto(flow_transformed, aes(x = !!sym(fluor_channels[1]), 
                                      y = !!sym(fluor_channels[2]))) +
        geom_hex(bins = 128) + 
        ggtitle("Before Compensation") +
        theme_minimal()
    
    p2 <- ggcyto(flow_compensated, aes(x = !!sym(fluor_channels[1]), 
                                      y = !!sym(fluor_channels[2]))) +
        geom_hex(bins = 128) + 
        ggtitle("After Compensation") +
        theme_minimal()
    
    gridExtra::grid.arrange(p1, p2, ncol = 2)
}

# Automated gating for flow cytometry analysis
library(openCyto)
library(flowDensity)

# Step 1: Remove debris and doublets
# FSC-A vs SSC-A for debris removal
debris_gate <- fsApply(flow_compensated, function(fr) {
    flowDensity::deGate(fr, channel = c("FSC-A", "SSC-A"), 
                       percentile = 0.95, upper = TRUE)
})

# Apply debris removal
cells_clean <- Subset(flow_compensated, debris_gate)
cat("📊 Cells after debris removal:", nrow(cells_clean), "\n")

# Step 2: Singlet selection (FSC-H vs FSC-A)
if("FSC-H" %in% colnames(cells_clean)) {
    singlet_gate <- fsApply(cells_clean, function(fr) {
        # Create singlet gate based on FSC-H vs FSC-A correlation
        flowDensity::deGate(fr, channel = c("FSC-A", "FSC-H"), 
                           percentile = 0.95, upper = TRUE)
    })
    
    singlets <- Subset(cells_clean, singlet_gate)
    cat("📊 Singlet cells:", nrow(singlets), "\n")
} else {
    singlets <- cells_clean
}

# Step 3: Live/Dead discrimination (if viability dye present)
viability_channels <- grep("Live|Dead|Viability|7AAD|PI", 
                          colnames(singlets), value = TRUE)

if(length(viability_channels) > 0) {
    live_gate <- fsApply(singlets, function(fr) {
        flowDensity::deGate(fr, channel = viability_channels[1], 
                           percentile = 0.95, upper = FALSE)
    })
    
    live_cells <- Subset(singlets, live_gate)
    cat("📊 Live cells:", nrow(live_cells), "\n")
} else {
    live_cells <- singlets
    cat("ℹ️ No viability dye detected\n")
}

# Step 4: Major population identification
# Example: CD3+ T cells vs CD19+ B cells
if(all(c("CD3", "CD19") %in% colnames(live_cells))) {
    # T cell gate (CD3+)
    t_cell_gate <- fsApply(live_cells, function(fr) {
        flowDensity::deGate(fr, channel = "CD3", percentile = 0.05)
    })
    
    # B cell gate (CD19+)  
    b_cell_gate <- fsApply(live_cells, function(fr) {
        flowDensity::deGate(fr, channel = "CD19", percentile = 0.05)
    })
    
    t_cells <- Subset(live_cells, t_cell_gate)
    b_cells <- Subset(live_cells, b_cell_gate)
    
    cat("🔬 Population Analysis Results:\n")
    cat("T cells (CD3+):", nrow(t_cells), 
        paste0("(", round(nrow(t_cells)/nrow(live_cells)*100, 1), "%)\n"))
    cat("B cells (CD19+):", nrow(b_cells), 
        paste0("(", round(nrow(b_cells)/nrow(live_cells)*100, 1), "%)\n"))
}

# Visualize gating strategy
if(length(fluor_channels) >= 2) {
    gating_plot <- ggcyto(live_cells, aes(x = !!sym(fluor_channels[1]), 
                                         y = !!sym(fluor_channels[2]))) +
        geom_hex(bins = 128) +
        ggtitle("Flow Cytometry Analysis: Final Gated Population") +
        theme_minimal() +
        labs(x = fluor_channels[1], y = fluor_channels[2])
    
    print(gating_plot)
}

# Statistical analysis for flow cytometry data
library(broom)
library(effsize)

# Create sample metadata (adjust based on your experimental design)
sample_data <- data.frame(
    sample_id = paste0("Sample_", 1:10),
    condition = rep(c("Control", "Treatment"), each = 5),
    replicate = rep(1:5, 2),
    stringsAsFactors = FALSE
)

# Calculate population frequencies for each sample
# This example assumes you have multiple FCS files
calculate_frequencies <- function(fcs_file, sample_info) {
    # Load and process each file
    flow_data <- read.FCS(fcs_file)
    # ... apply same preprocessing steps ...
    
    # Calculate frequencies
    total_cells <- nrow(flow_data)
    frequencies <- list(
        t_cells = nrow(t_cells) / total_cells * 100,
        b_cells = nrow(b_cells) / total_cells * 100,
        # Add other populations...
    )
    
    return(data.frame(
        sample_id = sample_info$sample_id,
        condition = sample_info$condition,
        frequencies
    ))
}

# Example with simulated data for demonstration
set.seed(123)
frequency_data <- data.frame(
    sample_id = sample_data$sample_id,
    condition = sample_data$condition,
    t_cells = c(rnorm(5, 45, 5), rnorm(5, 52, 5)),  # Control vs Treatment
    b_cells = c(rnorm(5, 15, 3), rnorm(5, 12, 3)),
    nk_cells = c(rnorm(5, 8, 2), rnorm(5, 11, 2))
)

# Statistical testing for each population
populations <- c("t_cells", "b_cells", "nk_cells")
statistical_results <- data.frame()

for(pop in populations) {
    # Perform t-test
    test_result <- t.test(frequency_data[[pop]] ~ frequency_data$condition)
    
    # Calculate effect size (Cohen's d)
    effect_size <- cohen.d(frequency_data[[pop]], frequency_data$condition)
    
    # Store results
    result_row <- data.frame(
        population = pop,
        p_value = test_result$p.value,
        mean_control = test_result$estimate[1],
        mean_treatment = test_result$estimate[2],
        cohens_d = effect_size$estimate,
        stringsAsFactors = FALSE
    )
    
    statistical_results <- rbind(statistical_results, result_row)
}

# Apply multiple testing correction
statistical_results$p_adjusted <- p.adjust(statistical_results$p_value, 
                                         method = "fdr")

# Identify significant changes
statistical_results$significant <- statistical_results$p_adjusted < 0.05
statistical_results$effect_magnitude <- ifelse(abs(statistical_results$cohens_d) > 0.8, "Large",
                                              ifelse(abs(statistical_results$cohens_d) > 0.5, "Medium", "Small"))

cat("📊 Statistical Analysis Results:\n")
print(statistical_results)

# Visualize results
library(ggplot2)
frequency_long <- reshape2::melt(frequency_data, 
                                id.vars = c("sample_id", "condition"),
                                variable.name = "population",
                                value.name = "frequency")

comparison_plot <- ggplot(frequency_long, 
                         aes(x = condition, y = frequency, fill = condition)) +
    geom_boxplot(alpha = 0.7) +
    geom_jitter(width = 0.2, alpha = 0.6) +
    facet_wrap(~population, scales = "free_y") +
    scale_fill_manual(values = c("Control" = "#3498db", "Treatment" = "#e74c3c")) +
    labs(title = "Flow Cytometry Analysis: Population Frequencies",
         subtitle = "Comparison between experimental conditions",
         x = "Condition", y = "Frequency (%)",
         fill = "Condition") +
    theme_minimal() +
    theme(strip.text = element_text(face = "bold"))

print(comparison_plot)