---

title: "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

---

# Project Description

This project aims to identify therapeutic targets for melanoma patients by following a 2-phase approach:

-   Phase I: Analyze scRNA-seq data to identify DE genes and infer enriched pathways\

-   Phase II: Compare responders vs. non-responders to identify therapeutic targets

## Dataset Description

This analysis integrates five complementary single-cell RNA sequencing datasets:

1. Dataset 1 (Treatment Response Dataset)

   - Contains treatment response data

   - Includes responder/non-responder classifications

   - Focus on immunotherapy outcomes

2. Dataset 2 (Tumor Cell Dataset)

   - Characterizes tumor cell heterogeneity

   - Includes cell type annotations

   - Focus on tumor cell populations

3. Dataset 3 (Tumor Microenvironment Dataset)

   - Maps the tumor microenvironment

   - Contains malignant/non-malignant classifications

   - Includes detailed immune cell typing

4. Dataset 4 (Control Melanocyte Dataset)

   - Normal melanocyte control data

   - Reference for non-malignant state

   - Baseline expression profiles

5. Dataset 5 (Control Immune Dataset)

   - Normal immune cell control data

   - Reference for immune cell states

   - Baseline immune profiles

## Strategy

Use `seurat` and `harmony` to complete the following tasks:

1.  Load/process/combine data\

2.  Clustering and biomarker determination\

3.  Comparisons: disease vs. control

# Initialize environment

Install required packages.

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

# Define packages to install

pkg.list = c('svDialogs', 'vroom', 'dplyr', 'DT', 'Seurat', 'harmony', 

             'patchwork', 'ggplot2', 'ggrepel', 'openxlsx', 'progress')

# Define packages not already installed

pkg.install <- pkg.list[!(pkg.list %in% installed.packages()[, 'Package'])]

# Install uninstalled packages

if (length(pkg.install) > 0) {

  install.packages(pkg.install)

}

```

Load installed packages.

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

# Load packages

library(svDialogs)    # for prompting user-input

library(vroom)        # for quickly reading data

library(dplyr)        # for data processing

library(DT)           # to display datatables

library(Seurat)       # for scRNA-seq analysis

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(openxlsx)     # to write data to excel

library(progress)     # to display progress bar

```

Define custom functions.

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

# qc_filter(): filter data based on QC for scRNA-seq data

qc_filter <- function(obj, feat.t=c(200, 2500), pct.mt.t=5, var.method='vst', 

                      feat.n=2000, qc.plot=T, top.n=10, title='') {

  ############################ FUNCTION DESCRIPTION ############################

  # feat.t = lower and upper limits on unique gene counts

  # pct.mt.t = threshold of level in mitochondrial contamination

  # var.method = method for selecting highly variable genes

  # feat.n = number of variable genes to select

  # qc.plot = boolean whether to generate plots to decide downstream analyses

  # title = string to use for plot title

  ############################## BEGIN FUNCTION ################################

  # determine percentage of mitochondrial contamination

  obj[['pct.mt']] <- PercentageFeatureSet(obj, pattern='^MT-')

  # filter + nomalize + scale data

  obj <- obj %>% 

    subset(., subset=(nFeature_RNA > feat.t[1]) & (nFeature_RNA < feat.t[2]) & 

             (pct.mt < pct.mt.t)) %>% NormalizeData(.) %>% 

    FindVariableFeatures(., selection.method=var.method) %>% 

    ScaleData(.) %>% RunPCA(., features=VariableFeatures(object=.))

  # generate follow-up QC plots (if prompted)

  if (qc.plot) { 

    p1 <- VariableFeaturePlot(obj) %>% 

      LabelPoints(plot=., points=head(VariableFeatures(obj), top.n), repel=T)

    p2 <- ElbowPlot(obj)

    plot(p1 + p2 + plot_annotation(title=title))

  }

  # return output

  return(obj)

}

# de_analyze(): conduct differential expression (DE) analysis

de_analyze <- function(m1, m2, alt='two.sided', paired=F, var.equal=F, 

                       adj.method='bonferroni', t=0.05, de.plot=F, title='') { 

  ############################ FUNCTION DESCRIPTION ############################

  # m1, m2 = expression matrices to compare

  # alt, paired, var.equal = arguments for t.test() function

  # adj.method = method for calculating adjusted p-value

  # t = threshold for significance 

  # de.plot = boolean whether to generate a volcano plot

  ############################## BEGIN FUNCTION ################################

  # make sure two matrices have same number of rows

  if (nrow(m1) != nrow(m2)) { 

    stop('Row length does not match between the provided matrices.')

  }

  # make sure gene names align between matrices

  if (!all(rownames(m1) == rownames(m2))) { 

    stop('Gene names do not align between provided matrices.')

    }

  # instantiate output variable

  results <- data.frame(gene=rownames(m1), 

                        t.stat=vector(mode='numeric', length=nrow(m1)), 

                        p.val=vector(mode='numeric', length=nrow(m1)))

  # conduct unpaired t-test with unequal variance for each gene

  pb <- progress_bar$new(

    format='  analyzing [:bar] :percent time left: :eta', total=nrow(m1))

  for (i in 1:nrow(m1)) { 

    pb$tick()

    x <- m1[i, ]; y <- m2[i, ]

    r <- t.test(x, y, alternative=alt, paired=paired, var.equal=var.equal)

    results$t.stat[i] <- r$statistic

    results$p.val[i] <- r$p.value

  }

  # determine adjusted p-values

  results$q.val <- p.adjust(results$p.val, method=adj.method)

  # add additional fields

  results <- results %>%

    mutate(Significance=case_when(q.val < t & t.stat > 0 ~ 'Up',

                                  q.val < t & t.stat < 0 ~ 'Down',

                                  T ~ 'NS')) %>% arrange(q.val)

  # generate volcano plot (if prompted)

  if (de.plot) { 

    p <- results %>% arrange(t.stat) %>% ggplot(data=., 

           aes(x=t.stat, y=-log10(q.val), col=Significance, label=gene)) +

      geom_point() + geom_text_repel() + theme_minimal() + 

      scale_color_manual(values=c('blue', 'black', 'red')) + ggtitle(title)

    plot(p)

  }

  # return output

  return(results)

}

```

Additional settings.

```{r settings}

# Adjust system settings

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

# Save plots? (default: F)

save.plots <- dlgInput('Save all outputs? (T/F)', F)$res

```

# Tasks

## 1. Load/process/combine data

**Description**: load all scRNA-seq files.

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

setwd("..")

# Load workspace (if it exists)

f <- list.files()

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

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

}

```

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

##################### ONLY RUN IF THERE'S NO WORKSPACE #########################

# # Load disease datasets

# obj <- list(); meta <- list()

# #   Dataset 1: Treatment response dataset

# meta$dataset1 <- read.delim(file='data/dataset1/metadata.txt')

# x <- vroom('data/dataset1/counts.txt', col_names=F)

# x <- x[, 1:nrow(meta$dataset1)] %>% as.matrix(.)

# y <- read.delim('data/dataset1/genes.txt', header=F)

# rownames(x) <- make.unique(y$V1)

# colnames(x) <- meta$dataset1$ID

# obj$dataset1 <- CreateSeuratObject(x, project='dataset1',

#                                     min.cells=3, min.features=200)

# obj$dataset1$Treatment <- as.factor(meta$dataset1$Treatment)

# obj$dataset1$Response <- as.factor(meta$dataset1$Response)

# rm(x, y); gc()

# 

# #   Dataset 2: Tumor cell dataset

# meta$dataset2 <- read.csv('data/dataset2/metadata.csv')

# x <- vroom('data/dataset2/expression_matrix.csv')

# y <- as.matrix(x[, -1]); rownames(y) <- x$...1

# obj$dataset2 <- CreateSeuratObject(y, project='dataset2',

#                                     min.cells=3, min.features=200)

# obj$dataset2$type <- as.factor(meta$dataset2$cell.types)

# 

# #   Dataset 3: Tumor microenvironment dataset

# x <- vroom('data/dataset3/expression_matrix.txt')

# y <- as.matrix(x[-c(1:3), -1]); rownames(y) <- x$Cell[-c(1:3)]

# obj$dataset3 <- CreateSeuratObject(y, project='dataset3',

#                                    min.cells=3, min.features=200)

# meta$dataset3 <- data.frame(Tumor=t(x[1, -1]),

#                             Malignant=case_when(x[2, -1] == 1 ~ 'No',

#                                                 x[2, -1] == 2 ~ 'Yes',

#                                                 x[2, -1] == 0 ~ 'Unresolved'),

#                             NMType=case_when(x[3, -1] == 1 ~ 'T',

#                                              x[3, -1] == 2 ~ 'B',

#                                              x[3, -1] == 3 ~ 'Macro',

#                                              x[3, -1] == 4 ~ 'Endo',

#                                              x[3, -1] == 5 ~ 'CAF',

#                                              x[3, -1] == 6 ~ 'NK',

#                                              x[3, -1] == 0 ~ 'NA'))

# obj$dataset3$Malignant <- meta$dataset3$Malignant

# obj$dataset3$NMType <- meta$dataset3$NMType

# 

# # Load control datasets

# #   Dataset 4: Control melanocyte dataset

# x <- vroom('data/dataset4/expression_matrix.csv')

# y <- as.matrix(x[, -1]); rownames(y) <- x$...1

# obj$dataset4 <- CreateSeuratObject(y, project='dataset4',

#                                     min.cells=3, min.features=200)

# 

# #   Dataset 5: Control immune dataset

# meta$dataset5 <- read.table('data/dataset5/metadata.tsv',

#                              sep='\t', header=T)

# x <- vroom('data/dataset5/matrix.tsv')

# y <- as.matrix(x[, -1]); rownames(y) <- x$Gene.Name

# obj$dataset5 <- CreateSeuratObject(y, project='dataset5',

#                                     min.cells=3, min.features=200)

# obj$dataset5$sample <- as.factor(meta$dataset5$sample)

# rm(f, x, y, meta); gc()

```

*Description*: Combine all datasets together into single Seurat object. Also apply `harmony` to remove clustering bias based on dataset source.

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

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

setwd(file_path)

# Combine all datasets

if (!exists('data.all')) { 

  data.all <- list()

}

if (!('raw' %in% names(data.all))) { 

  data.all$raw <- merge(obj$dataset1, obj$dataset2) %>%

    merge(., obj$dataset3) %>% merge(., obj$dataset4) %>% merge(., obj$dataset5)

}

# Add source information

if (!('Source' %in% names(data.all$raw@meta.data))) { 

  data.all$raw$Source <- c(rep(names(obj)[1], ncol(obj$dataset1)),

                           rep(names(obj)[2], ncol(obj$dataset2)),

                           rep(names(obj)[3], ncol(obj$dataset3)),

                           rep(names(obj)[4], ncol(obj$dataset4)),

                           rep(names(obj)[5], ncol(obj$dataset5)))

}

# Apply QC

if (!('proc' %in% names(data.all))) { 

  data.all$proc <- qc_filter(data.all$raw)

}

# Visualize integration results

p1 <- DimPlot(object=data.all$proc, reduction='pca', pt.size=0.1, 

              group.by='Source')

p2 <- VlnPlot(object=data.all$proc, features='PC_1', pt.size=0.1, 

              group.by='Source')

p1 + p2

if (save.plots) { 

  ggsave('plots/QC_no_harmony.pdf', width=10, height=5, units='in')

}

# Apply harmony (to remove clustering based on dataset source)

if (!('harmony' %in% names(data.all$proc@reductions))) { 

  data.all$proc <- data.all$proc %>% RunHarmony('Source', plot_convergence=T)

}

# Visualize integration results (after harmony)

p1 <- DimPlot(object=data.all$proc, reduction='harmony', pt.size=0.1, 

              group.by='Source')

p2 <- VlnPlot(object=data.all$proc, features='harmony_1', pt.size=0.1, 

              group.by='Source')

p1 + p2

if (save.plots) { 

  ggsave('plots/QC_w_harmony.pdf', width=10, height=5, units='in')

}

```

## 2. Clustering and biomarker determination

**Description**: cluster cells based on UMAP and determine biomarkers based on cluster assignment.

### Clustering

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

# Visualize UMAP plots (runtime: ~5 minutes)

if (!('umap' %in% names(data.all$proc@reductions))) { 

  data.all$proc <- data.all$proc %>% 

      RunUMAP(reduction='harmony', dims=1:20) %>% 

      FindNeighbors(reduction='harmony', dims=1:20) %>% 

      FindClusters(resolution=0.5) %>% identity()

}

#   by dataset source

p <- DimPlot(data.all$proc, reduction='umap', group.by='Source', pt.size=0.1, 

        split.by='Source') + ggtitle('UMAP split by dataset source'); plot(p)

#   by cluster (unlabeled)

p <- DimPlot(data.all$proc, reduction='umap', label=T, pt.size=0.1) + 

  ggtitle('UMAP of combined scRNA-seq data (unlabeled)'); plot(p)

if (save.plots) { 

  ggsave('plots/UMAP_unlabeled.pdf', width=10, height=5, units='in')

}

#   by cluster (labeled)

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)

p <- DimPlot(plot.data, reduction='umap', label=T, pt.size=0.1) + 

  ggtitle('UMAP of combined scRNA-seq data (labeled)'); plot(p)

if (save.plots) { 

  ggsave('plots/UMAP_labeled.pdf', width=10, height=5, units='in')

}

```

### Biomarkers

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

# Find all biomarkers based on clustering (runtime: ~30 minutes)

if (!('bm' %in% names(data.all))) { 

  data.all$bm <- FindAllMarkers(data.all$proc, min.pct=0.25, logfc.threshold=0.25)

}

# View table of top 3 biomarkers for each cluster

d <- data.all$bm %>% group_by(cluster) %>% slice_max(n=3, order_by=avg_log2FC)

datatable(d)

# Visualize ridge plots based on biomarkers of interest

ridge.feat <- dlg_list(sort(unique(d$gene)), multiple=T)$res

p <- RidgePlot(data.all$proc, features=ridge.feat, 

          ncol=ceiling(length(ridge.feat) / 2)); plot(p)

if (save.plots) { 

  ggsave('plots/biomarker_ridge_plots.pdf', width=10, height=10, units='in')

}

```

## 3. Comparisons: disease vs. control

**Description**: determine differentially expressed genes (DEGs) between disease vs. control groups.

### Case: tumor

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

# Instantiate DE variable

if (!exists('de')) { 

  de <- list()

}

if (!('data' %in% names(de))) { 

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

}

# Define datasets to compare

ix <- intersect(rownames(obj$dataset3), rownames(obj$dataset2)) %>%

  intersect(., rownames(obj$dataset4))

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

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

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

# DE analysis (runtime: ~15 minutes)

if (!('tumor' %in% names(de))) { 

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

}

# Visualize dot plot

dot.feat <- dlg_list(de$tumor$gene[de$tumor$q.val < 0.05], multiple=T)$res

p.data <- data.all$proc 

Idents(p.data) <- as.factor(case_when(jx1 ~ 'Disease', jx2 ~ 'Control', T ~ 'NA'))

p <- p.data %>% 

  DotPlot(object=., features=dot.feat, idents=c('Disease', 'Control')) + 

  ggtitle('Case: tumor')

plot(p)

if (save.plots) { 

  ggsave('plots/DE_tumor.pdf', width=10, height=5, units='in')

}

```

### Case: immune cells (bulk)

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

# Define datasets to compare

ix <- intersect(rownames(obj$dataset1), rownames(obj$dataset2)) %>% 

  intersect(., rownames(obj$dataset5))

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

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

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

# DE analysis (runtime: ~20 minutes)

if (!('immune' %in% names(de))) { 

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

}

# Visualize dot plot

dot.feat <- dlg_list(de$immune$gene[de$immune$q.val < 0.05], multiple=T)$res

p.data <- data.all$proc 

Idents(p.data) <- as.factor(case_when(jx1 ~ 'Disease', jx2 ~ 'Control', T ~ 'NA'))

p <- p.data %>% 

  DotPlot(object=., features=dot.feat, idents=c('Disease', 'Control')) + 

  ggtitle('Case: immune cells')

plot(p)

if (save.plots) { 

  ggsave('plots/DE_immune_bulk.pdf', width=10, height=5, units='in')

}

```

### Case: immune cells (cluster-based)

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

# Define datasets to compare (runtime: ~5 minutes)

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

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

  # Define row + column indices

  ix <- intersect(rownames(obj$dataset1), rownames(obj$dataset2)) %>% 

    intersect(., rownames(obj$dataset5))

  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

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

  # DE analysis

  if (!(immune.cells[i] %in% names(de))) { 

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

  }

  # Visualize dot plot

  x <- de[[immune.cells[[i]]]]

  dot.feat <- x$gene[x$q.val < 0.05][1:5]

  p.data <- data.all$proc 

  Idents(p.data) <- as.factor(case_when(jx1 ~ 'Disease', jx2 ~ 'Control', T ~ 'NA'))

  p <- p.data %>% 

    DotPlot(object=., features=dot.feat, idents=c('Disease', 'Control')) + 

    ggtitle(paste0('Case: ', immune.cells[i]))

  plot(p)

  if (save.plots) { 

    ggsave(paste0('plots/DE_', immune.cells[i], '.pdf'), 

           width=10, height=5, units='in')

  }

}

```

# Case: responder vs. non-responder

ix <- rownames(obj$dataset1)

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

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

# Save outputs and working space (if prompted)

```

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

# Define key glycolytic genes

glycolysis_genes <- c(

  "HK1", "HK2", "GPI", "PFKL", "PFKM", "PFKP",  # Glucose uptake and early glycolysis

  "ALDOA", "ALDOB", "ALDOC",                     # Aldolase isoforms

  "TPI1",                                        # Triose phosphate isomerase

  "GAPDH",                                       # Glyceraldehyde phosphate dehydrogenase

  "PGK1",                                        # Phosphoglycerate kinase

  "PGAM1",                                       # Phosphoglycerate mutase

  "ENO1", "ENO2",                               # Enolase isoforms

  "PKM", "PKLR",                                # Pyruvate kinase isoforms

  "LDHA", "LDHB",                               # Lactate dehydrogenase

  "SLC2A1", "SLC2A3"                            # Glucose transporters (GLUT1, GLUT3)

)

# Filter for genes present in the dataset

glycolysis_genes_present <- intersect(glycolysis_genes, rownames(de$data))

# Create anonymous gene IDs

anon_genes <- paste0("gene_", seq_along(glycolysis_genes_present))

names(anon_genes) <- glycolysis_genes_present

# Create mapping dataframe for reference

gene_mapping <- data.frame(

  Original_Gene = glycolysis_genes_present,

  Anonymous_ID = anon_genes[glycolysis_genes_present]

)

# Write mapping to file if saving is enabled

if (save.plots) {

  write.csv(gene_mapping, 'results/glycolysis_gene_mapping.csv', row.names = FALSE)

}

# Create dot plot for glycolytic genes

p.data <- data.all$proc

Idents(p.data) <- as.factor(case_when(jx1 ~ 'Disease', jx2 ~ 'Control', T ~ 'NA'))

# Generate dot plot with anonymous IDs and custom colors

p_glycolysis <- DotPlot(p.data, 

                        features = glycolysis_genes_present,

                        idents = c('Disease', 'Control')) +

  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +

  xlab('Genes') +

  ylab('Condition') +

  scale_color_gradient2(low = "navy", mid = "gray90", high = "red", 

                       midpoint = 0, name = "Average\nExpression") +

  scale_x_discrete(labels = anon_genes[glycolysis_genes_present])

# Display plot

print(p_glycolysis)

# Save plot if requested

if (save.plots) {

  ggsave('plots/DE_tumor_glycolysis.pdf', 

         plot = p_glycolysis, 

         width = 12, 

         height = 6, 

         units = 'in')

}

# Calculate statistics for glycolytic genes

glycolysis_stats <- de$tumor[de$tumor$gene %in% glycolysis_genes_present, ]

glycolysis_stats <- glycolysis_stats[order(glycolysis_stats$q.val), ]

# Add anonymous IDs to statistics

glycolysis_stats$anonymous_id <- anon_genes[glycolysis_stats$gene]

# Display statistics with anonymous IDs

print("Differential Expression Statistics for De-identified Glycolytic Genes:")

glycolysis_stats$gene <- glycolysis_stats$anonymous_id

print(glycolysis_stats)

# Count significantly different genes

sig_glyco_genes <- sum(glycolysis_stats$q.val < 0.05)

print(sprintf("Number of significantly different glycolytic genes (q < 0.05): %d", sig_glyco_genes))

```

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

# Get expression data for glycolytic genes

expr_data <- GetAssayData(p.data, slot = "data")[glycolysis_genes_present,]

# Filter out genes with all zeros

nonzero_genes <- rownames(expr_data)[rowSums(expr_data > 0) > 0]

expr_data <- expr_data[nonzero_genes,]

# Update anonymous gene IDs for non-zero genes

anon_genes <- anon_genes[nonzero_genes]

target_genes <- c("g10", "g21", "g03")

selected_genes <- names(anon_genes)[anon_genes %in% target_genes]

# Create a data frame for plotting and store statistical results

plot_data <- data.frame()

stat_results <- data.frame(

  Gene = character(),

  P_value = numeric(),

  T_statistic = numeric(),

  Mean_diff = numeric(),

  stringsAsFactors = FALSE

)

for (gene in selected_genes) {

  # Get expression values for current gene

  expr_values <- expr_data[gene,]

  # Split values by condition before centering

  disease_vals <- expr_values[colnames(expr_data) %in% colnames(de$data)[jx1]]

  control_vals <- expr_values[colnames(expr_data) %in% colnames(de$data)[jx2]]

  # Perform t-test on raw values

  t_test_result <- t.test(disease_vals, control_vals)

  # Center the expression values by subtracting the overall mean

  overall_mean <- mean(expr_values)

  disease_vals_centered <- disease_vals - overall_mean

  control_vals_centered <- control_vals - overall_mean

  # Store statistical results

  stat_results <- rbind(stat_results, data.frame(

    Gene = anon_genes[gene],

    P_value = t_test_result$p.value,

    T_statistic = t_test_result$statistic,

    Mean_diff = mean(disease_vals) - mean(control_vals),

    Mean_centered_diff = mean(disease_vals_centered) - mean(control_vals_centered),

    stringsAsFactors = FALSE

  ))

  # Create temporary data frames for each condition with centered values

  temp_df <- data.frame(

    Expression = c(disease_vals_centered, control_vals_centered),

    Gene = anon_genes[gene],

    Condition = c(rep("Disease", length(disease_vals_centered)), 

                 rep("Control", length(control_vals_centered)))

  )

  # Append to main data frame

  plot_data <- rbind(plot_data, temp_df)

}

# Add FDR correction

stat_results$FDR <- p.adjust(stat_results$P_value, method = "BH")

# Sort results by p-value

stat_results <- stat_results[order(stat_results$P_value), ]

# Create boxplot with swarm plot

p_boxplot <- ggplot(plot_data, aes(x = Gene, y = Expression, fill = Condition)) +

  # Add swarm plot layer first (so it appears behind)

  geom_jitter(aes(color = Condition), position = position_jitterdodge(jitter.width = 0.2), 

              size = 0.3, alpha = 0.3) +

  # Add boxplot layer second (so it appears in front)

  geom_boxplot(outlier.shape = NA, alpha = 0.4, width = 0.7) +

  # Colors for boxplot fill and point colors

  scale_fill_manual(values = c("Disease" = "#E41A1C", "Control" = "#377EB8")) +

  scale_color_manual(values = c("Disease" = "#E41A1C", "Control" = "#377EB8")) +

  # Set y-axis limits

  ylim(-2, max(plot_data$Expression)) +

  # Theme customization

  theme_minimal() +

  theme(

    axis.text.x = element_text(angle = 45, hjust = 1),

    panel.grid.major.y = element_blank(),  # Remove y-axis major grid lines

    panel.grid.minor.y = element_blank(),  # Remove y-axis minor grid lines

    panel.grid.major.x = element_line(color = "gray90"),  # Keep x-axis major grid lines

    panel.grid.minor.x = element_blank(),  # Remove x-axis minor grid lines

    legend.position = c(0.95, 0.95),  # Position legend inside top-right

    legend.justification = c(1, 1),    # Anchor point for legend

    legend.box.just = "right",

    legend.margin = margin(6, 6, 6, 6),

    legend.background = element_rect(fill = "white", color = NA),

    legend.box.background = element_rect(color = "gray90")

  ) +

  # Labels

  xlab("Genes") +

  ylab("Centered Expression") +

  ggtitle("Gene Expression Relative to Mean") +

  # Add horizontal line at y=0

  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50", alpha = 0.5)

# Display plot

print(p_boxplot)

# Save plot if requested

if (save.plots) {

  ggsave('plots/DE_tumor_glycolysis_boxplots.pdf', 

         plot = p_boxplot, 

         width = 12, 

         height = 6, 

         units = 'in')

}

# Calculate summary statistics for centered data

summary_stats <- plot_data %>%

  group_by(Gene, Condition) %>%

  summarise(

    Mean = mean(Expression),

    Median = median(Expression),

    SD = sd(Expression),

    Q1 = quantile(Expression, 0.25),

    Q3 = quantile(Expression, 0.75),

    n = n(),

    .groups = 'drop'

  )

# Display summary statistics

print("Summary Statistics for Centered Expression by Condition:")

print(summary_stats)

# Display statistical test results

print("\nStatistical Test Results (Disease vs Control):")

print(stat_results)

# Print number of genes removed due to zero expression

n_removed <- length(glycolysis_genes_present) - length(nonzero_genes)

if (n_removed > 0) {

  print(sprintf("\nRemoved %d genes with zero expression", n_removed))

  print("Removed genes:")

  print(setdiff(glycolysis_genes_present, nonzero_genes))

}

# Print number of significantly different genes

sig_genes <- sum(stat_results$FDR < 0.05)

if (sig_genes > 0) {

  print(sprintf("\nNumber of significantly different genes (FDR < 0.05): %d", sig_genes))

  print("Significant genes:")

  print(stat_results[stat_results$FDR < 0.05, ])

}

```