---

title: "DE Analysis of scRNA-seq data"

date: "Last modified: `r format(Sys.time(), '%B %d, %Y')`"

tags: [scRNA-seq, seurat, melanoma, immunotherapy, PBMCs] 

output:

  html_document:

    theme: flatly

    highlight: zenburn

    toc: true

    number_sections: false

    toc_depth: 2

    toc_float:

      collapsed: false

      smooth_scroll: true

    code_folding: hide

    self_contained: yes

---

# Notebook Description

This notebook compiles the code and outputs for differential expression (DE) analysis. Specifically, two methods are covered:  

  * unpaired Student's t-test (custom function) 

  * [MAST](https://doi.org/10.1186/s13059-015-0844-5) 

Of note, pre-processed data stored in the `scRNAseq_analysis.RData` file is used for the analyses detailed below.  

# Initialize environment

```{r load_packages, results="hide", message=F, warning=F, error=F}

# Install required packages if not already installed

if (!require("BiocManager", quietly = TRUE)) {

  install.packages("BiocManager")

}

if (!require("MAST", quietly = TRUE)) {

  BiocManager::install("MAST")

}

if (!require("presto", quietly = TRUE)) {

  install.packages("presto")

}

# Memory management settings

options(future.globals.maxSize = 16000 * 1024^2)  # Set maximum global size to 8GB

# Load core packages first

library(Matrix)       # for sparse matrix operations

library(Seurat)      # for scRNA-seq analysis

library(MAST)        # for scRNAseq DEG analysis

# Load other packages

library(svDialogs)    # for prompting user-input

library(vroom)        # for quickly reading data

library(dplyr)        # for data processing

library(DT)           # to display datatables

library(harmony)      # to integration scRNA-seq data

library(patchwork)    # for combining plots

library(ggplot2)      # for data visualization

library(ggrepel)      # to use geom_pont_repel()

library(ggvenn)       # to visualize venn diagrams

library(openxlsx)     # to write data to excel

library(progress)     # to display progress bar

library(rio)          # to load all worksheets in a workbook

library(pheatmap)     # to visualize heatmaps

library(ggfortify)    # to visualize PCA plots

library(presto)       # for faster Wilcoxon test implementation

library(gridExtra)    # for arranging multiple plots

# Set up parallel processing

library(future)

plan(multisession, workers = 2)  # Adjust based on available memory

options(future.globals.maxSize = 16000 * 1024^2)

```

System settings.  

```{r settings}

# Adjust system settings

Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 2)

# Enhanced memory management functions

clear_memory <- function(deep_clean = FALSE) {

  # Standard cleanup

  gc(verbose = FALSE)

  if(exists(".Random.seed")) rm(.Random.seed)

  # Deep cleanup if requested

  if(deep_clean) {

    # Clear package caches

    if(exists(".__C__")) rm(list=".__C__", envir=.GlobalEnv)

    # Clear hidden objects

    rm(list = ls(all.names = TRUE, envir = .GlobalEnv), envir = .GlobalEnv)

    # Force garbage collection multiple times

    replicate(3, gc(verbose = FALSE))

  }

  invisible(gc())

}

# Function to get data efficiently with batching support

get_assay_data <- function(object, data_type = "data", batch_size = 1000) {

  tryCatch({

    # Get data using appropriate method

    if (packageVersion("Seurat") >= "5.0.0") {

      data <- GetAssayData(object = object, layer = data_type)

    } else {

      data <- GetAssayData(object = object, slot = data_type)

    }

    # Convert to sparse matrix if dense

    if (!is(data, "dgCMatrix")) {

      data <- as(data, "dgCMatrix")

    }

    return(data)

  }, error = function(e) {

    stop(sprintf("Failed to get assay data: %s", e$message))

  })

}

# Function to process data in batches

process_in_batches <- function(data, batch_size = 1000, FUN, ...) {

  n_cells <- ncol(data)

  n_batches <- ceiling(n_cells / batch_size)

  results <- vector("list", n_batches)

  pb <- progress_bar$new(total = n_batches)

  for(i in 1:n_batches) {

    start_idx <- (i-1) * batch_size + 1

    end_idx <- min(i * batch_size, n_cells)

    # Process batch

    batch_data <- data[, start_idx:end_idx]

    results[[i]] <- FUN(batch_data, ...)

    # Clean up

    rm(batch_data)

    clear_memory()

    pb$tick()

  }

  # Combine results

  combined_results <- do.call(rbind, results)

  rm(results)

  clear_memory()

  return(combined_results)

}

# Add gene anonymization function after the settings chunk

anonymize_genes <- function(gene_list, prefix = "GENE") {

  # Create mapping of real genes to anonymous IDs

  anon_ids <- paste0(prefix, "_", sprintf("%04d", seq_along(gene_list)))

  names(anon_ids) <- gene_list

  return(anon_ids)

}

```

Load pre-processed data. 

```{r load_data, warning=F, message=F}

file_path <- "/Users/scampit/Projects/torchstack/cancer-biomarker-discovery"

setwd(file_path)

# Start with clean environment

#clear_memory(deep_clean = TRUE)

# Load data from saved workspace

f <- list.files()

if (any(endsWith(f, '.RData'))) {

  load(f[endsWith(f, '.RData')][1])

} else {

  stop("No .RData file found. Please ensure scRNAseq_wksp.RData is in the working directory.")

}

# Initialize DE lists

de <- list()

# Get data from Seurat object efficiently

tryCatch({

  message("Loading expression data from Seurat object...")

  de$data <- get_assay_data(data.all$proc)

  # Verify data loaded correctly

  if (is.null(de$data) || !is(de$data, "dgCMatrix")) {

    stop("Failed to load expression data properly")

  }

  message(sprintf("Loaded expression matrix with %d genes and %d cells", 

                 nrow(de$data), ncol(de$data)))

}, error = function(e) {

  stop(sprintf("Error loading data: %s", e$message))

})

#clear_memory()

# Load or initialize results

de$ttest <- if(file.exists('results/tables/DE_ttest.xlsx')) {

  import_list('results/tables/DE_ttest.xlsx')

} else {

  list()

}

de$mast <- if(file.exists('results/tables/DE_mast.xlsx')) {

  import_list('results/tables/DE_mast.xlsx')

} else {

  list()

}

#clear_memory(deep_clean = TRUE)

# Verify data structure

str(de, max.level = 1)

```

```{r summary_stats, warning=F, message=F}

# Assign row names to MAST results

f <- names(de$ttest)

for (i in 1:length(f)) { 

  x <- de$ttest[[f[i]]]; y <- de$mast[[f[i]]]

  g <- intersect(x$gene[x$q.val < 0.05], y$gene[y$p_val_adj < 0.05]) 

  cat(sprintf('%s\tDEGs (t-test): %d\tDEGs (MAST): %d\tDEGs (overlap): %d\n', 

              f[i], sum(x$q.val < 0.05), sum(y$p_val_adj < 0.05), length(g)))

}

```

# DE analyses

**Description**: determine differentially expressed genes (DEGs) between different comparisons.  

```{r DE_initialize, warning=F, message=T, eval=!('data' %in% de)}

# Instantiate DE variable + relevant fields

de <- list(); de$ttest <- list(); de$mast <- list()

# Define data to work with - update to use layer

tryCatch({

  de$data <- GetAssayData(object=data.all$proc, layer='data')

}, error = function(e) {

  # Fallback for older Seurat versions

  message("Attempting fallback for older Seurat version...")

  de$data <- GetAssayData(object=data.all$proc, slot='data')

})

# Add validation

if (is.null(de$data) || ncol(de$data) == 0 || nrow(de$data) == 0) {

  stop("Failed to retrieve data from Seurat object")

}

```

## Case: disease vs. control (tumor)

Comparison between benign vs. tumor skin cells. 

Estimated runtime: ~15 minutes. 

```{r DE_tumor, warning=F, message=T, eval=!('tumor' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Define cell groups to compare

ix <- intersect(rownames(obj$GSE72056), rownames(obj$GSE115978)) %>%

  intersect(., rownames(obj$GSE151091))

jx1 <- data.all$proc$Source == names(obj)[3] | data.all$proc$Malignant %in% 'Yes'

jx2 <- data.all$proc$Source == names(obj)[4]

# Subsample cells if too many

max_cells_per_group <- 1000

if(sum(jx1) > max_cells_per_group) {

  set.seed(42)

  cells_g1 <- sample(which(jx1), max_cells_per_group)

  jx1[setdiff(which(jx1), cells_g1)] <- FALSE

}

if(sum(jx2) > max_cells_per_group) {

  set.seed(42)

  cells_g2 <- sample(which(jx2), max_cells_per_group)

  jx2[setdiff(which(jx2), cells_g2)] <- FALSE

}

# DE analysis (t-test) with memory management

x <- de$data[ix, jx1]

y <- de$data[ix, jx2]

#clear_memory()

de$ttest$tumor <- de_analyze(x, y) %>% na.omit

rm(x, y)

#clear_memory(deep_clean = TRUE)

# DE analysis (MAST)

message("Creating Seurat object...")

combined_data <- cbind(de$data[ix, jx1], de$data[ix, jx2])

if (!is(combined_data, "dgCMatrix")) {

  combined_data <- as(combined_data, "dgCMatrix")

}

x <- CreateSeuratObject(counts = combined_data)

rm(combined_data)

#clear_memory()

# Process data in steps with memory management

message("Normalizing data...")

x <- NormalizeData(x, verbose = FALSE)

#clear_memory()

message("Finding variable features...")

x <- FindVariableFeatures(x, verbose = FALSE)

#clear_memory()

message("Scaling data...")

x <- ScaleData(x, verbose = FALSE)

#clear_memory()

message("Running PCA...")

x <- RunPCA(x, verbose = FALSE)

#clear_memory()

# Set identities

cells_g1 <- 1:sum(jx1)

cells_g2 <- (sum(jx1) + 1):ncol(x)

Idents(object = x, cells = cells_g1) <- 'Disease'

Idents(object = x, cells = cells_g2) <- 'Control'

rm(cells_g1, cells_g2)

#clear_memory()

# Run MAST analysis with error handling and memory management

message("Running differential expression analysis...")

tryCatch({

  de$mast$tumor <- FindMarkers(

    object = x,

    ident.1 = 'Disease',

    ident.2 = 'Control',

    test.use = 'MAST',

    verbose = TRUE,

    max.cells.per.ident = 1000,

    min.pct = 0.1,        # Filter low-expressed genes

    logfc.threshold = 0.25 # Filter low effect size

  )

}, error = function(e) {

  message("Error in MAST analysis: ", e$message)

  message("Attempting alternative analysis method...")

  de$mast$tumor <- FindMarkers(

    object = x,

    ident.1 = 'Disease',

    ident.2 = 'Control',

    test.use = 'wilcox',

    verbose = TRUE,

    max.cells.per.ident = 1000,

    min.pct = 0.1,

    logfc.threshold = 0.25

  )

})

# Clean up

rm(x)

#clear_memory(deep_clean = TRUE)

# End timer + log time elapsed

t.end <- Sys.time()

message("Analysis completed in: ", difftime(t.end, t.start, units = "mins"), " minutes")

```

## Case: disease vs. control (immune cells, bulk)

Bulk comparison between healthy vs. diseased immune cells. 

Estimated runtime: ~20 minutes. 

```{r DE_immune_bulk,warning=F, message=T, eval=!('immune' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Define cell groups to compare

ix <- intersect(rownames(obj$GSE120575), rownames(obj$GSE115978)) %>% 

  intersect(., rownames(obj$GSE183212))

jx1 <- data.all$proc$Source == names(obj)[1] | data.all$proc$Malignant %in% 'No'

jx2 <- data.all$proc$Source == names(obj)[5]

# DE analysis (t-test)

x <- de$data[ix, jx1]; y <- de$data[ix, jx2]

de$ttest$immune <- de_analyze(x, y) %>% na.omit()

# DE analysis (MAST)

# Create Seurat object with proper assay

x <- CreateSeuratObject(counts = cbind(de$data[ix, jx1], de$data[ix, jx2]))

# Normalize the data

x <- NormalizeData(x)

# Find variable features

x <- FindVariableFeatures(x)

# Scale the data

x <- ScaleData(x)

# Run PCA

x <- RunPCA(x)

# Set identities

Idents(object=x, cells=1:sum(jx1)) <- 'Tumor'

Idents(object=x, cells=sum(jx1)+1:ncol(x)) <- 'Benign'

# Run MAST analysis

de$mast$immune <- FindMarkers(object=x, 

                             ident.1='Tumor', 

                             ident.2='Benign', 

                             test.use='MAST',

                             verbose = TRUE)

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

## Case: disease vs. control (immune cells, cluster-based)

Immune cell-specific comparisons between healthy vs. diseased cells. 

Estimated runtime: ~15 minutes. 

```{r DE_immune_cluster, warning=F, message=T}

# Start timer

t.start <- Sys.time()

# Define labeled scRNAseq data

lab <- c('Tumor cell', 'Technical error 1', 'Cytotoxic CD8 T cell 1', 'B cell', 

         'Melanocytes 1', 'Macrophage', 'Cytotoxic CD8 T cell 2', 'Technical error 2', 

         'Memory T cell', 'Melanocytes 2', 'Dysfunctional CD8 T cell', 

         'Cytotoxic CD8 T cell 3', '? 1', 'Melanocytes 3', 'Activated cell', 

         'Exhausted T cell', 'Cytotoxic T cell', '? 2', '? 3', 'Melanocytes 4', '? 4')

names(lab) <- levels(data.all$proc)

plot.data <- data.all$proc %>% RenameIdents(., lab)

# Make sure we have the expression data

if (!exists("de") || is.null(de$data) || !is(de$data, "dgCMatrix")) {

  message("Reloading expression data...")

  de <- list()

  de$data <- GetAssayData(object = data.all$proc, layer = "data")

  if (!is(de$data, "dgCMatrix")) {

    de$data <- as(de$data, "dgCMatrix")

  }

}

# Initialize results lists if they don't exist

if (!"ttest" %in% names(de)) de$ttest <- list()

if (!"mast" %in% names(de)) de$mast <- list()

# Iterate through each immune cell group

immune.cells <- c('B', 'CD8', 'Macrophage', 'Memory', 'Activated')

for (i in 1:length(immune.cells)) { 

  message(sprintf("\nProcessing %s cells...", immune.cells[i]))

  # Define row + column indices

  ix <- intersect(rownames(obj$GSE120575), rownames(obj$GSE115978)) %>% 

    intersect(., rownames(obj$GSE183212))

  c <- grepl(immune.cells[i], Idents(plot.data))

  jx1 <- (plot.data$Source == names(obj)[1] | 

            data.all$proc$Malignant %in% 'No') & c

  jx2 <- data.all$proc$Source == names(obj)[5] & c

  # Print number of genes + cells

  message(sprintf('No. of genes: %d\tNo. of G1: %d\tNo. of G2: %d', 

                length(ix), sum(jx1), sum(jx2)))

  # Skip if not enough cells

  if (sum(jx1) < 3 || sum(jx2) < 3) {

    message(sprintf("Skipping %s - insufficient cells", immune.cells[i]))

    next

  }

  # DE analysis (t-test)

  message("Running t-test analysis...")

  x <- de$data[ix, jx1]

  y <- de$data[ix, jx2]

  de$ttest[[immune.cells[i]]] <- de_analyze(x, y) %>% na.omit()

  rm(x, y)

  clear_memory()

  # DE analysis (MAST)

  message("Running MAST analysis...")

  tryCatch({

    # Create Seurat object

    combined_data <- cbind(de$data[ix, jx1], de$data[ix, jx2])

    if (!is(combined_data, "dgCMatrix")) {

      combined_data <- as(combined_data, "dgCMatrix")

    }

    x <- CreateSeuratObject(counts = combined_data)

    rm(combined_data)

    clear_memory()

    # Process data

    x <- x %>%

      NormalizeData(verbose = FALSE) %>%

      FindVariableFeatures(verbose = FALSE) %>%

      ScaleData(verbose = FALSE)

    # Set identities

    Idents(object = x, cells = 1:sum(jx1)) <- 'Tumor'

    Idents(object = x, cells = (sum(jx1) + 1):ncol(x)) <- 'Benign'

    # Run MAST

    de$mast[[immune.cells[i]]] <- FindMarkers(

      object = x,

      ident.1 = 'Tumor',

      ident.2 = 'Benign',

      test.use = 'MAST',

      verbose = TRUE,

      max.cells.per.ident = 1000

    )

    rm(x)

    clear_memory()

  }, error = function(e) {

    message(sprintf("Error in MAST analysis for %s: %s", immune.cells[i], e$message))

    message("Attempting Wilcoxon test instead...")

    # Try Wilcoxon as fallback

    de$mast[[immune.cells[i]]] <- FindMarkers(

      object = x,

      ident.1 = 'Tumor',

      ident.2 = 'Benign',

      test.use = 'wilcox',

      verbose = TRUE,

      max.cells.per.ident = 1000

    )

  })

  message(sprintf("Completed analysis for %s", immune.cells[i]))

}

# End timer + log time elapsed

t.end <- Sys.time()

message("\nTotal analysis time: ", difftime(t.end, t.start, units = "mins"), " minutes")

# Print summary of results

message("\nSummary of analyses:")

for (cell_type in immune.cells) {

  if (cell_type %in% names(de$ttest)) {

    n_ttest <- nrow(de$ttest[[cell_type]])

    n_mast <- if(cell_type %in% names(de$mast)) nrow(de$mast[[cell_type]]) else 0

    message(sprintf("%s: %d t-test DEGs, %d MAST DEGs", cell_type, n_ttest, n_mast))

  }

}

```

## Case: responder vs. non-responder (GSE120575)

Comparison between responder vs. non-responder cells to immunotherapy. 

Estimated runtime: ~1 hour.  

```{r DE_response, warning=F, message=T, eval=!('response' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Make sure we have the expression data

if (!exists("de") || is.null(de$data) || !is(de$data, "dgCMatrix")) {

  message("Reloading expression data...")

  de <- list()

  de$data <- GetAssayData(object = data.all$proc, layer = "data")

  if (!is(de$data, "dgCMatrix")) {

    de$data <- as(de$data, "dgCMatrix")

  }

}

# Define cell groups to compare

ix <- rownames(obj$GSE120575)

jx1 <- data.all$proc$Response %in% 'Responder'

jx2 <- data.all$proc$Response %in% 'Non-responder'

# Print cell counts

message(sprintf("Number of responder cells: %d", sum(jx1)))

message(sprintf("Number of non-responder cells: %d", sum(jx2)))

# Subsample cells if too many

max_cells_per_group <- 1000

if(sum(jx1) > max_cells_per_group) {

  set.seed(42)

  cells_g1 <- sample(which(jx1), max_cells_per_group)

  jx1[setdiff(which(jx1), cells_g1)] <- FALSE

}

if(sum(jx2) > max_cells_per_group) {

  set.seed(42)

  cells_g2 <- sample(which(jx2), max_cells_per_group)

  jx2[setdiff(which(jx2), cells_g2)] <- FALSE

}

# DE analysis (t-test)

message("Running t-test analysis...")

x <- de$data[ix, jx1]

y <- de$data[ix, jx2]

de$ttest$response <- de_analyze(x, y) %>% na.omit

rm(x, y)

#clear_memory()

# DE analysis (MAST)

message("Running MAST analysis...")

tryCatch({

  # Create Seurat object

  message("Creating Seurat object...")

  combined_data <- cbind(de$data[ix, jx1], de$data[ix, jx2])

  if (!is(combined_data, "dgCMatrix")) {

    combined_data <- as(combined_data, "dgCMatrix")

  }

  x <- CreateSeuratObject(counts = combined_data)

  rm(combined_data)

  clear_memory()

  # Process data

  message("Processing data...")

  x <- x %>%

    NormalizeData(verbose = FALSE) %>%

    FindVariableFeatures(verbose = FALSE) %>%

    ScaleData(verbose = FALSE) %>%

    RunPCA(verbose = FALSE)

  # Set identities

  message("Setting identities...")

  cells_g1 <- 1:sum(jx1)

  cells_g2 <- (sum(jx1) + 1):ncol(x)

  Idents(object = x, cells = cells_g1) <- 'Responder'

  Idents(object = x, cells = cells_g2) <- 'Non-responder'

  # Verify identities

  if (!all(levels(Idents(x)) %in% c('Responder', 'Non-responder'))) {

    stop("Failed to set identities correctly")

  }

  # Run MAST

  message("Running differential expression analysis...")

  de$mast$response <- FindMarkers(

    object = x,

    ident.1 = 'Responder',

    ident.2 = 'Non-responder',

    test.use = 'MAST',

    verbose = TRUE,

    max.cells.per.ident = 1000,

    min.pct = 0.1,

    logfc.threshold = 0.25

  )

  # Clean up

  rm(x, cells_g1, cells_g2)

  clear_memory()

}, error = function(e) {

  message("Error in MAST analysis: ", e$message)

  message("Attempting Wilcoxon test instead...")

  # Try Wilcoxon as fallback

  de$mast$response <- FindMarkers(

    object = x,

    ident.1 = 'Responder',

    ident.2 = 'Non-responder',

    test.use = 'wilcox',

    verbose = TRUE,

    max.cells.per.ident = 1000,

    min.pct = 0.1,

    logfc.threshold = 0.25

  )

})

# End timer + log time elapsed

t.end <- Sys.time()

message("\nAnalysis completed in: ", difftime(t.end, t.start, units = "mins"), " minutes")

# Print summary

if (!is.null(de$mast$response)) {

  message("\nFound ", sum(de$mast$response$p_val_adj < 0.05), 

          " significant DEGs (FDR < 0.05)")

}

```

# Visualizations

**Description**: visual inspection of DE analysis results. 

## T-test vs. MAST DEGs

Visual comparison of DEGs determined with t-test vs. MAST (venn diagrams).

```{r DE_venn, warning=F, message=T}

# Compare DEGs b/t the two methods

f <- names(de$mast); p <- list()

for (i in 1:length(f)) { 

  x <- de$ttest[[f[i]]]; y <- de$mast[[f[i]]]

  g <- intersect(x$gene[x$q.val < 0.05], row.names(y)[y$p_val_adj < 0.05])

  # g <- intersect(x$gene[x$q.val < 0.05], y$gene[y$p_val_adj < 0.05]) 

  cat(sprintf('No. of intersecting DEGs for %s: %d\n', f[i], length(g)))

  ggvenn(list('t-test'=x$gene[x$q.val < 0.05], 

              'MAST'=row.names(y)[y$p_val_adj < 0.05])) + ggtitle(f[i])

  ggsave(paste0('plots/DE_venn_', f[i], '.pdf'), width=5, height=5, units='in')

}

```

## Responder vs. non-responder DEGs

Clustergram (or heatmap) of DEGs identified for responder vs. non-responder comparison. 

```{r DE_response_heat, warning=F, message=T}

# Modify DE_response_heat chunk

tryCatch({

  # Get significant genes

  if (!exists("de") || is.null(de$mast$response)) {

    stop("MAST results not found. Please run the DE analysis first.")

  }

  sig_genes <- rownames(de$mast$response)[de$mast$response$p_val_adj < 0.05]

  if (length(sig_genes) == 0) {

    stop("No significant genes found (FDR < 0.05)")

  }

  message(sprintf("Found %d significant genes", length(sig_genes)))

  # Take top N genes by absolute log fold change

  n <- min(50, length(sig_genes))

  g <- sig_genes[order(abs(de$mast$response$avg_log2FC[match(sig_genes, rownames(de$mast$response))]), 

                      decreasing = TRUE)][1:n]

  # Create anonymous gene IDs

  anon_genes <- anonymize_genes(g)

  # Get all cells

  jx1 <- which(data.all$proc$Response %in% 'Responder')

  jx2 <- which(data.all$proc$Response %in% 'Non-responder')

  # Get expression data

  x <- de$data[g, c(jx1, jx2)]

  if (!is(x, "dgCMatrix")) {

    x <- as(x, "dgCMatrix")

  }

  x <- as.matrix(x)

  # Replace gene names with anonymous IDs

  rownames(x) <- anon_genes[rownames(x)]

  # Calculate mean expression for each group

  x_resp <- rowMeans(x[, 1:length(jx1), drop=FALSE])

  x_nonresp <- rowMeans(x[, (length(jx1)+1):ncol(x), drop=FALSE])

  # Combine into matrix

  x_means <- cbind(x_resp, x_nonresp)

  colnames(x_means) <- c("Responder\n(mean)", "Non-responder\n(mean)")

  # Calculate standard errors for error bars

  x_resp_se <- apply(x[, 1:length(jx1), drop=FALSE], 1, function(x) sd(x)/sqrt(length(x)))

  x_nonresp_se <- apply(x[, (length(jx1)+1):ncol(x), drop=FALSE], 1, function(x) sd(x)/sqrt(length(x)))

  # Add SE information to rownames

  rownames(x_means) <- sprintf("%s\n(SE: %.2f, %.2f)", 

                              rownames(x_means), 

                              x_resp_se, 

                              x_nonresp_se)

  # Scale data

  x_scaled <- t(scale(t(x_means)))

  # Create annotation

  col.annot <- data.frame(

    Group = c('Responder', 'Non-responder'),

    row.names = colnames(x_means)

  )

  # Color palette

  col.pal <- colorRampPalette(colors = c('navy', 'white', 'red'))(100)

  # Create heatmap

  message("Generating heatmap...")

  p <- pheatmap(

    x_scaled,

    cluster_rows = TRUE,

    cluster_cols = FALSE,

    scale = "none",  # already scaled above

    color = col.pal,

    annotation_col = col.annot,

    show_rownames = TRUE,

    show_colnames = TRUE,

    main = sprintf('Top %d De-identified DEGs (Mean Expression)', nrow(x_scaled)),

    fontsize_row = 8,

    annotation_colors = list(

      Group = c('Responder' = '#E41A1C', 'Non-responder' = '#377EB8')

    )

  )

  # Save plot

  message("Saving plot...")

  pdf('plots/DE_heatmap_response_collapsed.pdf', width = 11, height = 8)

  print(p)

  dev.off()

  # Save gene ID mapping

  write.csv(data.frame(

    Anonymous_ID = names(anon_genes),

    Original_Gene = anon_genes,

    Mean_Responder = x_resp,

    SE_Responder = x_resp_se,

    Mean_NonResponder = x_nonresp,

    SE_NonResponder = x_nonresp_se

  ), './response_gene_mapping_with_stats.csv', row.names = FALSE)

  message("Heatmap generated successfully")

  message("Gene mapping with statistics saved to results/response_gene_mapping_with_stats.csv")

}, error = function(e) {

  message(sprintf("Error generating heatmap: %s", e$message))

})

```

# Add new chunk after DE_response_heat

```{r DE_response_heat_individual, warning=F, message=T}

# Show individual cell heatmap alongside collapsed view

tryCatch({

  # Get significant genes

  if (!exists("de") || is.null(de$mast$response)) {

    stop("MAST results not found. Please run the DE analysis first.")

  }

  sig_genes <- rownames(de$mast$response)[de$mast$response$p_val_adj < 0.05]

  if (length(sig_genes) == 0) {

    stop("No significant genes found (FDR < 0.05)")

  }

  message(sprintf("Found %d significant genes", length(sig_genes)))

  # Take top N genes by absolute log fold change

  n <- min(50, length(sig_genes))

  g <- sig_genes[order(abs(de$mast$response$avg_log2FC[match(sig_genes, rownames(de$mast$response))]), 

                      decreasing = TRUE)][1:n]

  # Create anonymous gene IDs

  anon_genes <- anonymize_genes(g)

  # Get cells (subsample to make visualization manageable)

  jx1 <- which(data.all$proc$Response %in% 'Responder')

  jx2 <- which(data.all$proc$Response %in% 'Non-responder')

  # Subsample cells

  cells_per_group <- min(25, min(length(jx1), length(jx2)))

  set.seed(42)

  sampled_jx1 <- sample(jx1, cells_per_group)

  sampled_jx2 <- sample(jx2, cells_per_group)

  # Get expression data for individual cells

  x_individual <- de$data[g, c(sampled_jx1, sampled_jx2)]

  if (!is(x_individual, "dgCMatrix")) {

    x_individual <- as(x_individual, "dgCMatrix")

  }

  x_individual <- as.matrix(x_individual)

  # Clean individual cell data

  message("Cleaning individual cell data...")

  # Replace Inf/-Inf with NA

  x_individual[is.infinite(x_individual)] <- NA

  # Remove rows (genes) with any NA values

  complete_rows <- complete.cases(x_individual)

  if (!all(complete_rows)) {

    message(sprintf("Removing %d genes with NA values", sum(!complete_rows)))

    x_individual <- x_individual[complete_rows, ]

    g <- g[complete_rows]

    anon_genes <- anon_genes[complete_rows]

  }

  # Check for zero variance genes

  gene_var <- apply(x_individual, 1, var)

  if (any(gene_var == 0)) {

    message(sprintf("Removing %d genes with zero variance", sum(gene_var == 0)))

    x_individual <- x_individual[gene_var > 0, ]

    g <- g[gene_var > 0]

    anon_genes <- anon_genes[gene_var > 0]

  }

  # Replace gene names with anonymous IDs

  rownames(x_individual) <- anon_genes[rownames(x_individual)]

  # Scale individual cell data with validation

  message("Scaling individual cell data...")

  x_individual_scaled <- t(scale(t(x_individual)))

  # Check for NA/Inf values after scaling

  if (any(is.na(x_individual_scaled)) || any(is.infinite(x_individual_scaled))) {

    message("Found NA/Inf values after scaling, capping extreme values at ±10")

    x_individual_scaled[x_individual_scaled > 10] <- 10

    x_individual_scaled[x_individual_scaled < -10] <- -10

    x_individual_scaled[is.na(x_individual_scaled)] <- 0

  }

  # Create annotation for individual cells

  col.annot.individual <- data.frame(

    Response = rep(c('Responder', 'Non-responder'), each = cells_per_group),

    row.names = colnames(x_individual)

  )

  # Get mean expression data (for side-by-side comparison)

  x_resp <- rowMeans(de$data[g, jx1, drop=FALSE])

  x_nonresp <- rowMeans(de$data[g, jx2, drop=FALSE])

  x_means <- cbind(x_resp, x_nonresp)

  colnames(x_means) <- c("Responder\n(mean)", "Non-responder\n(mean)")

  # Calculate standard errors

  x_resp_se <- apply(de$data[g, jx1, drop=FALSE], 1, function(x) sd(x)/sqrt(length(x)))

  x_nonresp_se <- apply(de$data[g, jx2, drop=FALSE], 1, function(x) sd(x)/sqrt(length(x)))

  # Clean mean data

  message("Cleaning mean data...")

  x_means[is.infinite(x_means)] <- NA

  complete_rows_means <- complete.cases(x_means)

  if (!all(complete_rows_means)) {

    message(sprintf("Removing %d genes with NA values from means", sum(!complete_rows_means)))

    x_means <- x_means[complete_rows_means, ]

    x_resp_se <- x_resp_se[complete_rows_means]

    x_nonresp_se <- x_nonresp_se[complete_rows_means]

  }

  # Add SE information to rownames for mean data

  rownames(x_means) <- sprintf("%s\n(SE: %.2f, %.2f)", 

                              anon_genes[rownames(x_means)], 

                              x_resp_se, 

                              x_nonresp_se)

  # Scale mean data with validation

  message("Scaling mean data...")

  x_means_scaled <- t(scale(t(x_means)))

  if (any(is.na(x_means_scaled)) || any(is.infinite(x_means_scaled))) {

    message("Found NA/Inf values after scaling means, capping extreme values at ±10")

    x_means_scaled[x_means_scaled > 10] <- 10

    x_means_scaled[x_means_scaled < -10] <- -10

    x_means_scaled[is.na(x_means_scaled)] <- 0

  }

  # Create annotation for means

  col.annot.means <- data.frame(

    Group = c('Responder', 'Non-responder'),

    row.names = colnames(x_means)

  )

  # Color palette

  col.pal <- colorRampPalette(colors = c('navy', 'white', 'red'))(100)

  message("Generating plots...")

  # Create individual cells heatmap

  p1 <- pheatmap(

    x_individual_scaled,

    cluster_rows = TRUE,

    cluster_cols = FALSE,

    scale = "none",

    color = col.pal,

    annotation_col = col.annot.individual,

    show_rownames = TRUE,

    show_colnames = FALSE,

    gaps_col = cells_per_group,

    main = sprintf('Individual Cells\n(n=%d per group)', cells_per_group),

    fontsize_row = 8,

    annotation_colors = list(

      Response = c('Responder' = '#E41A1C', 'Non-responder' = '#377EB8')

    ),

    legend = FALSE,  # Hide colorbar

    treeheight_row = 0,  # Hide row dendrogram

    cellwidth = 10,  # Force square cells

    cellheight = 10  # Force square cells

  )

  # Create mean expression heatmap

  p2 <- pheatmap(

    x_means_scaled,

    cluster_rows = TRUE,

    cluster_cols = FALSE,

    scale = "none",

    color = col.pal,

    annotation_col = col.annot.means,

    show_rownames = TRUE,

    show_colnames = TRUE,

    main = 'Group Means\n(all cells)',

    fontsize_row = 8,

    annotation_colors = list(

      Group = c('Responder' = '#E41A1C', 'Non-responder' = '#377EB8')

    ),

    legend = FALSE,  # Hide colorbar

    treeheight_row = 0,  # Hide row dendrogram

    cellwidth = 20,  # Force square cells (larger since fewer columns)

    cellheight = 20,  # Force square cells (larger since fewer columns)

    labels_row = sub("\n.*", "", rownames(x_means))  # Only show gene IDs without SE info

  )

  # Save plots to PDF

  pdf('plots/DE_heatmap_response_comparison.pdf', width = 15, height = 8)

  grid::grid.newpage()

  grid::grid.draw(gridExtra::grid.arrange(p1$gtable, p2$gtable, ncol=2))

  dev.off()

  # Display plots in notebook

  grid::grid.newpage()

  grid::grid.draw(gridExtra::grid.arrange(p1$gtable, p2$gtable, ncol=2))

  # Save detailed statistics

  write.csv(data.frame(

    Anonymous_ID = names(anon_genes),

    Original_Gene = anon_genes,

    Mean_Responder = x_resp,

    SE_Responder = x_resp_se,

    Mean_NonResponder = x_nonresp,

    SE_NonResponder = x_nonresp_se,

    N_Responder = length(jx1),

    N_NonResponder = length(jx2)

  ), 'results/response_gene_mapping_detailed.csv', row.names = FALSE)

  message("Side-by-side heatmaps generated successfully")

  message("Detailed gene statistics saved to results/response_gene_mapping_detailed.csv")

}, error = function(e) {

  message(sprintf("Error generating heatmaps: %s", e$message))

  message("\nDebugging information:")

  message("1. Check if expression data contains extreme values")

  message("2. Verify gene filtering and scaling steps")

  message("3. Ensure all matrices have proper dimensions")

  message("4. Look for NA/Inf values in the data")

})

```

PCA of DEGs identified for responder vs. non-responder comparison. 

```{r DE_response_PCA, warning=F, message=T}

# Modify DE_response_PCA chunk

tryCatch({

  # Validate MAST results exist

  if (!exists("de") || is.null(de$mast$response)) {

    stop("MAST results not found. Please run the DE analysis first.")

  }

  # Get significant genes

  sig_genes <- rownames(de$mast$response)[de$mast$response$p_val_adj < 0.05]

  if (length(sig_genes) == 0) {

    stop("No significant genes found (FDR < 0.05)")

  }

  message(sprintf("Found %d significant genes", length(sig_genes)))

  # Get responder/non-responder cells

  jx <- data.all$proc$Response %in% c('Responder', 'Non-responder')

  if (sum(jx) == 0) {

    stop("No cells found matching response criteria")

  }

  message(sprintf("Found %d cells (%d Responders, %d Non-responders)", 

                 sum(jx),

                 sum(data.all$proc$Response[jx] == "Responder"),

                 sum(data.all$proc$Response[jx] == "Non-responder")))

  # Subsample cells if too many

  max_cells_per_group <- 1000

  responder_idx <- which(data.all$proc$Response[jx] == "Responder")

  nonresponder_idx <- which(data.all$proc$Response[jx] == "Non-responder")

  if (length(responder_idx) > max_cells_per_group) {

    set.seed(42)

    responder_idx <- sample(responder_idx, max_cells_per_group)

  }

  if (length(nonresponder_idx) > max_cells_per_group) {

    set.seed(42)

    nonresponder_idx <- sample(nonresponder_idx, max_cells_per_group)

  }

  selected_cells <- c(responder_idx, nonresponder_idx)

  # Get expression data

  x <- de$data[sig_genes, jx][, selected_cells]

  if (!is(x, "dgCMatrix")) {

    x <- as(x, "dgCMatrix")

  }

  x <- as.matrix(x)

  # Create anonymous gene IDs

  anon_genes <- anonymize_genes(rownames(x))

  rownames(x) <- anon_genes[rownames(x)]

  # Handle zero variance genes

  gene_var <- apply(x, 1, var)

  if (any(gene_var == 0)) {

    message(sprintf("Removing %d genes with zero variance", sum(gene_var == 0)))

    x <- x[gene_var > 0, ]

  }

  # Clean data

  x[is.infinite(x)] <- NA

  x[is.na(x)] <- 0

  # Transpose and convert to data frame

  x_t <- t(x)

  # Create metadata

  m <- data.frame(

    Response = data.all$proc$Response[jx][selected_cells],

    row.names = rownames(x_t)

  )

  # Apply PCA with validation

  if (ncol(x_t) < 2) {

    stop("Need at least 2 features for PCA")

  }

  if (nrow(x_t) < 3) {

    stop("Need at least 3 samples for PCA")

  }

  message("Running PCA...")

  response_pca <- prcomp(x_t, scale. = TRUE)

  # Calculate variance explained

  var_explained <- summary(response_pca)$importance[2, 1:2] * 100

  # Calculate group centroids

  centroids <- aggregate(cbind(PC1, PC2) ~ Response, 

                        data = cbind(as.data.frame(response_pca$x[,1:2]), m), 

                        FUN = mean)

  # Create enhanced PCA plot

  p <- ggplot(data = as.data.frame(response_pca$x[,1:2]), 

              aes(x = PC1, y = PC2, color = m$Response)) +

    # Add points with reduced alpha for transparency

    geom_point(alpha = 0.3, size = 2) +

    # Add confidence ellipses

    stat_ellipse(level = 0.95, size = 1) +

    # Add centroids

    geom_point(data = centroids, aes(x = PC1, y = PC2), 

               size = 5, shape = 18) +

    # Add centroid labels

    geom_text(data = centroids, 

              aes(x = PC1, y = PC2, label = Response),

              vjust = -1.5, size = 4, color = "black") +