Datasets
Models
Downloads: 1
749 lines (624 with data), 27.5 kB
---
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, ])
}
```