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") +
# Customize appearance
ggtitle('De-identified DEG PCA with Group Trends') +
xlab(sprintf("PC1 (%.1f%%)", var_explained[1])) +
ylab(sprintf("PC2 (%.1f%%)", var_explained[2])) +
theme_minimal() +
scale_color_manual(values = c('Responder' = '#E41A1C', 'Non-responder' = '#377EB8')) +
theme(
legend.position = "bottom",
plot.title = element_text(hjust = 0.5, size = 14),
axis.title = element_text(size = 12),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
)
# Save plot
message("Saving plot...")
ggsave('plots/DE_PCA_response_collapsed.pdf', plot = p, width = 8, height = 8, units = 'in')
# Display plot
print(p)
message("PCA visualization completed successfully")
}, error = function(e) {
message(sprintf("Error in PCA analysis: %s", e$message))
message("Please ensure:")
message("1. DE analysis has been run successfully")
message("2. There are significant DEGs to analyze")
message("3. There are enough cells in each response group")
message("4. Expression data is valid and contains sufficient variation")
})
# 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))
})
# 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...")
# Set up PDF device for both plots
pdf('plots/DE_heatmap_response_comparison.pdf', width = 15, height = 8)
# Create layout for side-by-side plots
layout(matrix(1:2, ncol = 2))
# Plot 1: Individual cells
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
)
# Plot 2: Mean expression
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
)
dev.off()
# 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")
})
# 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)
# 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")
})