From b7a1efcdb25b8b763eec27ab7d0582106e61100f Mon Sep 17 00:00:00 2001 From: hsmaan Date: Wed, 27 Sep 2023 16:13:05 -0400 Subject: [PATCH] readding CIDER pbmc files --- ...l_Cider_Supp_Analysis_Figures_Stat_Tests.R | 1146 +++++++++++++++++ .../R/21_PBMC_perturbation_umap_plots.R | 206 +++ workflow/configs/config_control_cider.json | 22 + workflow/configs/config_umap.json | 50 + workflow/scripts/R/install_packages.R | 19 + .../scripts/python/cluster_classification.py | 135 ++ .../python/cluster_knn_classification.py | 135 ++ workflow/scripts/python/umap_plots.py | 129 ++ workflow/src/R/cider_integrate.R | 53 + workflow/src/python/utils/cider_integrate.py | 104 ++ workflow/src/python/utils/umap.py | 122 ++ 11 files changed, 2121 insertions(+) create mode 100644 workflow/analysis/R/17_Iniq_Contol_Cider_Supp_Analysis_Figures_Stat_Tests.R create mode 100644 workflow/analysis/R/21_PBMC_perturbation_umap_plots.R create mode 100644 workflow/configs/config_control_cider.json create mode 100644 workflow/configs/config_umap.json create mode 100644 workflow/scripts/R/install_packages.R create mode 100644 workflow/scripts/python/cluster_classification.py create mode 100644 workflow/scripts/python/cluster_knn_classification.py create mode 100644 workflow/scripts/python/umap_plots.py create mode 100644 workflow/src/R/cider_integrate.R create mode 100644 workflow/src/python/utils/cider_integrate.py create mode 100644 workflow/src/python/utils/umap.py diff --git a/workflow/analysis/R/17_Iniq_Contol_Cider_Supp_Analysis_Figures_Stat_Tests.R b/workflow/analysis/R/17_Iniq_Contol_Cider_Supp_Analysis_Figures_Stat_Tests.R new file mode 100644 index 0000000..7eb6ad5 --- /dev/null +++ b/workflow/analysis/R/17_Iniq_Contol_Cider_Supp_Analysis_Figures_Stat_Tests.R @@ -0,0 +1,1146 @@ +library(data.table) +library(tidyverse) +library(reshape2) +library(ggplot2) +library(ggthemes) +library(ggExtra) +library(ggpubr) +library(dotwhisker) +library(Seurat) +library(SeuratDisk) +library(ComplexHeatmap) +library(circlize) +library(RColorBrewer) +library(Cairo) +library(networkD3) + +# Helper functions +`%ni%` <- Negate(`%in%`) + +# Kevin's palette for plotting many catagoricals +kev_palette <- c( + "dodgerblue2", "#E31A1C", + "green4", + "#6A3D9A", + "#FF7F00", + "black", "gold1", + "skyblue2", "#FB9A99", + "palegreen2", + "#CAB2D6", + "#FDBF6F", + "gray70", "khaki2", + "maroon", "orchid1", "deeppink1", "blue1", "steelblue4", + "darkturquoise", "green1", "yellow4", "yellow3", + "darkorange4", "brown" +) + +# Change to results dir for control data w/ cider implemented +setwd("../../../results/control_w_cider/") + +##### Analysis of PBMC 2 batch balanced data - baseline ##### + +# Load in and concatenate imbalance summary files +setwd("imbalance_summaries/") +imba_files <- list.files() +imba_files <- grep( + "pbmc_2_batch_base_balanced", + imba_files, + value = TRUE +) +imba_loaded <- lapply(imba_files, fread) +imba_concat <- Reduce(rbind, imba_loaded) +gc() + +# Load in and concatenate the clustering summary results +setwd("../clustering_summaries/") +clus_files <- list.files() +clus_files <- grep( + "pbmc_2_batch_base_balanced", + clus_files, + value = TRUE +) +clus_loaded <- lapply(clus_files, fread) +clus_concat <- Reduce(rbind, clus_loaded) +clus_concat <- clus_concat[clus_concat$Method != "liger"] +gc() + +# Load in and concatenate clustering concordance summaries +setwd("../clustering_concord_summaries/") +clus_concord_files <- list.files() +clus_concord_files <- grep( + "pbmc_2_batch_base_balanced", + clus_concord_files, + value = TRUE +) +clus_concord_loaded <- lapply(clus_concord_files, fread) +clus_concord_concat <- Reduce(rbind, clus_concord_loaded) +clus_concord_concat <- clus_concord_concat[ + clus_concord_concat$`Method 1` != "liger" +] +clus_concord_concat <- clus_concord_concat[ + clus_concord_concat$`Method 2` != "liger" +] +gc() + +# Load in and concatenate cluster classification results +setwd("../clustering_classification/") +clus_class_files <- list.files() +clus_class_files <- grep( + "pbmc_2_batch_base_balanced", + clus_class_files, + value = TRUE +) +clus_class_loaded <- lapply(clus_class_files, fread) +clus_class_concat <- Reduce(rbind, clus_class_loaded) +clus_class_concat <- clus_class_concat[ + clus_class_concat$Method != "liger" +] +gc() + +# Load in and concatenate KNN-based cluster classification results +setwd("../knn_clustering_classification/") +knn_clus_class_files <- list.files() +knn_clus_class_files <- grep( + "pbmc_2_batch_base_balanced", + knn_clus_class_files, + value = TRUE +) +knn_clus_class_loaded <- lapply(knn_clus_class_files, fread) +knn_clus_class_concat <- Reduce(rbind, knn_clus_class_loaded) +knn_clus_class_concat <- knn_clus_class_concat[ + knn_clus_class_concat$Method != "liger" +] + +# Change to top level dir +setwd("../../../") + +# Create the cider output directories (figures/results) if it doesn't exist +if (!dir.exists("outs/control_w_cider/figures")) { + dir.create("outs/control_w_cider/figures", recursive = TRUE) +} +if (!dir.exists("outs/control_w_cider/results")) { + dir.create("outs/control_w_cider/results", recursive = TRUE) +} + +# Begin by plotting the ARI scores for each method + +# Merge imbalance and clustering summary results +imba_clus_merged <- merge( + clus_concat, + imba_concat, + by = c( + "Number of batches downsampled", + "Number of celltypes downsampled", + "Proportion downsampled", + "Replicate" + ) +) + +# Format celltype names +imba_clus_merged$`Downsampled celltypes` <- plyr::mapvalues( + imba_clus_merged$`Downsampled celltypes`, + from = c( + "Monocyte_CD14", + "Monocyte_FCGR3A", + "CD4 T cell", + "CD8 T cell" + ), + to = c( + "CD14+ Monocyte", + "FCGR3A+ Monocyte", + "CD4+ T cell", + "CD8+ T cell" + ) +) + +# Indicate which samples are controls and which are real runs +imba_clus_merged$type <- ifelse( + imba_clus_merged$`Number of batches downsampled` == 0, + "Control", + ifelse( + imba_clus_merged$`Proportion downsampled` == 0, + "Ablated", + "Downsampled" + ) +) + +# Get median celltype ARI based on each method and whether or not +# it's a control, downsampling, or ablation, and by celltype +median_celltype_ari_results <- imba_clus_merged %>% + group_by(Method, type, `Downsampled celltypes`) %>% + summarize( + `Median celltype ARI` = median(`Celltype ARI Imbalanced`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +# Melt and format for ComplexHeatMap plotting +median_celltype_ari_results_vals_long <- reshape2::dcast( + median_celltype_ari_results, + formula = type + `Downsampled celltypes` ~ `Method`, + value.var = "Median celltype ARI" +) +median_celltype_ari_results_vals_long$type <- factor( + median_celltype_ari_results_vals_long$type, + levels = c("Control", "Downsampled", "Ablated") +) +median_celltype_ari_results_vals_long <- median_celltype_ari_results_vals_long[ + order( + median_celltype_ari_results_vals_long$type, + median_celltype_ari_results_vals_long$`Downsampled celltypes` + ), +] +rownames(median_celltype_ari_results_vals_long) <- c( + paste0("R_", seq(1, nrow(median_celltype_ari_results_vals_long))) +) +colnames(median_celltype_ari_results_vals_long)[1] <- c( + "Type" +) +median_celltype_ari_long_vals_only <- median_celltype_ari_results_vals_long[ + ,-c(1,2) +] +median_celltype_ari_long_vals_only_scaled <- scale( + median_celltype_ari_long_vals_only, + center = TRUE, + scale = TRUE +) +median_celltype_ari_long_type <- median_celltype_ari_results_vals_long[ + ,1, drop = FALSE +] +median_celltype_ari_long_celltype <- median_celltype_ari_results_vals_long[ + ,2, drop = FALSE +] + +# Plot the three heatmaps together for median celltype ARI post integration +dark_2_cols = palette.colors(n = 8, "Dark2") +col_type = c( + "Control" = "forestgreen", + "Downsampled" = "darkorchid3", + "Ablated" = "firebrick2" +) +col_celltype = c( + "B cell" = dark_2_cols[1], + "CD14+ Monocyte" = dark_2_cols[2], + "CD4+ T cell" = dark_2_cols[3], + "CD8+ T cell" = dark_2_cols[4], + "FCGR3A+ Monocyte" = dark_2_cols[5], + "NK cell" = dark_2_cols[6], + "None" = "black" +) + +ht1 = Heatmap( + as.matrix(median_celltype_ari_long_vals_only_scaled), + name = "Scaled median \ncell-type ARI", + width = unit(5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE +) +ht2 = Heatmap( + as.matrix(median_celltype_ari_long_type), + name = "Type", + col = col_type, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +ht3 = Heatmap( + as.matrix(median_celltype_ari_long_celltype), + name = "Affected cell-type", + col = col_celltype, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +celltype_ari_hm <- ht1 + ht2 + ht3 +CairoPDF( + "outs/control_w_cider/figures/17_celltype_ari_ds_effects_heatmap_no_liger.pdf", + width = 8, + height = 6 +) +draw( + celltype_ari_hm, + column_title = "Integration method", + column_title_side = "bottom", + column_title_gp = gpar(fontsize = 14, fontface = "bold") +) +dev.off() + +# Perform the exact same analysis/heatmap as above, but for the balanced +# celltype ARI + +median_balanced_celltype_ari_results <- imba_clus_merged %>% + group_by(Method, type, `Downsampled celltypes`) %>% + summarize( + `Median balanced celltype ARI` = median(`Celltype ARI Balanced`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +# Melt and format for ComplexHeatMap plotting +median_balanced_celltype_ari_results_vals_long <- reshape2::dcast( + median_balanced_celltype_ari_results, + formula = type + `Downsampled celltypes` ~ `Method`, + value.var = "Median balanced celltype ARI" +) +median_balanced_celltype_ari_results_vals_long$type <- factor( + median_balanced_celltype_ari_results_vals_long$type, + levels = c("Control", "Downsampled", "Ablated") +) +median_balanced_celltype_ari_results_vals_long <- + median_balanced_celltype_ari_results_vals_long[ + order( + median_balanced_celltype_ari_results_vals_long$type, + median_balanced_celltype_ari_results_vals_long$`Downsampled celltypes` + ), + ] +rownames(median_balanced_celltype_ari_results_vals_long) <- c( + paste0("R_", seq(1, nrow(median_balanced_celltype_ari_results_vals_long))) + ) +colnames(median_balanced_celltype_ari_results_vals_long)[1] <- c( + "Type" +) +median_balanced_celltype_ari_long_vals_only <- + median_balanced_celltype_ari_results_vals_long[ + ,-c(1,2) + ] +median_balanced_celltype_ari_long_vals_only_scaled <- scale( + median_balanced_celltype_ari_long_vals_only, + center = TRUE, + scale = TRUE +) +median_balanced_celltype_ari_long_type <- + median_balanced_celltype_ari_results_vals_long[ + ,1, drop = FALSE + ] +median_balanced_celltype_ari_long_celltype <- + median_balanced_celltype_ari_results_vals_long[ + ,2, drop = FALSE + ] + +# Plot the three heatmaps together for median celltype ARI post integration +dark_2_cols = palette.colors(n = 8, "Dark2") +col_type = c( + "Control" = "forestgreen", + "Downsampled" = "darkorchid3", + "Ablated" = "firebrick2" +) +col_celltype = c( + "B cell" = dark_2_cols[1], + "CD14+ Monocyte" = dark_2_cols[2], + "CD4+ T cell" = dark_2_cols[3], + "CD8+ T cell" = dark_2_cols[4], + "FCGR3A+ Monocyte" = dark_2_cols[5], + "NK cell" = dark_2_cols[6], + "None" = "black" +) + +ht1 = Heatmap( + as.matrix(median_balanced_celltype_ari_long_vals_only_scaled), + name = "Scaled median \ncell-type balanced ARI", + width = unit(5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE +) +ht2 = Heatmap( + as.matrix(median_balanced_celltype_ari_long_type), + name = "Type", + col = col_type, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +ht3 = Heatmap( + as.matrix(median_balanced_celltype_ari_long_celltype), + name = "Affected cell-type", + col = col_celltype, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +celltype_ari_hm <- ht1 + ht2 + ht3 +CairoPDF( + "outs/control_w_cider/figures/17_celltype_balanced_ari_ds_effects_heatmap_no_liger.pdf", + width = 8, + height = 6 +) +draw( + celltype_ari_hm, + column_title = "Integration method", + column_title_side = "bottom", + column_title_gp = gpar(fontsize = 14, fontface = "bold") +) +dev.off() + +# Perform the exact same analysis/heatmap as above, but now for Batch ARI + +# Get median batch ARI based on each method and whether or not +# it's a control, downsampling, or ablation, and by celltype +median_batch_ari_results <- imba_clus_merged %>% + group_by(Method, type, `Downsampled celltypes`) %>% + summarize( + `Median batch ARI` = median(`Batch ARI`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +# Melt and format for ComplexHeatMap plotting +median_batch_ari_results_vals_long <- reshape2::dcast( + median_batch_ari_results, + formula = type + `Downsampled celltypes` ~ `Method`, + value.var = "Median batch ARI" +) +median_batch_ari_results_vals_long$type <- factor( + median_batch_ari_results_vals_long$type, + levels = c("Control", "Downsampled", "Ablated") +) +median_batch_ari_results_vals_long <- median_batch_ari_results_vals_long[ + order( + median_batch_ari_results_vals_long$type, + median_batch_ari_results_vals_long$`Downsampled celltypes` + ), +] +rownames(median_batch_ari_results_vals_long) <- c( + paste0("R_", seq(1, nrow(median_batch_ari_results_vals_long))) +) +colnames(median_batch_ari_results_vals_long)[1] <- c( + "Type" +) +median_batch_ari_long_vals_only <- median_batch_ari_results_vals_long[ + ,-c(1,2) +] +median_batch_ari_long_vals_only_scaled <- scale( + median_batch_ari_long_vals_only, + center = TRUE, + scale = TRUE +) +median_batch_ari_long_type <- median_batch_ari_results_vals_long[ + ,1, drop = FALSE +] +median_batch_ari_long_celltype <- median_batch_ari_results_vals_long[ + ,2, drop = FALSE +] + +# Plot the three heatmaps together for median batch ARI post integration +dark_2_cols = palette.colors(n = 8, "Dark2") +col_type = c( + "Control" = "forestgreen", + "Downsampled" = "darkorchid3", + "Ablated" = "firebrick2" +) +col_celltype = c( + "B cell" = dark_2_cols[1], + "CD14+ Monocyte" = dark_2_cols[2], + "CD4+ T cell" = dark_2_cols[3], + "CD8+ T cell" = dark_2_cols[4], + "FCGR3A+ Monocyte" = dark_2_cols[5], + "NK cell" = dark_2_cols[6], + "None" = "black" +) + +ht1 = Heatmap( + as.matrix(median_batch_ari_long_vals_only_scaled), + name = "Scaled median \nbatch ARI", + width = unit(5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE +) +ht2 = Heatmap( + as.matrix(median_batch_ari_long_type), + name = "Type", + col = col_type, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +ht3 = Heatmap( + as.matrix(median_batch_ari_long_celltype), + name = "Affected cell-type", + col = col_celltype, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +batch_ari_hm <- ht1 + ht2 + ht3 +CairoPDF( + "outs/control_w_cider/figures/17_batch_ari_ds_effects_heatmap_no_liger.pdf", + width = 8, + height = 6 +) +draw( + batch_ari_hm, + column_title = "Integration method", + column_title_side = "bottom", + column_title_gp = gpar(fontsize = 14, fontface = "bold") +) +dev.off() + +# Do the same analysis above for imbalanced cell-type ARI, but this time don't +# scale the ARI values - keep the absolute differences between methods +imba_clus_merged <- merge( + clus_concat, + imba_concat, + by = c( + "Number of batches downsampled", + "Number of celltypes downsampled", + "Proportion downsampled", + "Replicate" + ) +) + +imba_clus_merged$`Downsampled celltypes` <- plyr::mapvalues( + imba_clus_merged$`Downsampled celltypes`, + from = c( + "Monocyte_CD14", + "Monocyte_FCGR3A", + "CD4 T cell", + "CD8 T cell" + ), + to = c( + "CD14+ Monocyte", + "FCGR3A+ Monocyte", + "CD4+ T cell", + "CD8+ T cell" + ) +) + +imba_clus_merged$type <- ifelse( + imba_clus_merged$`Number of batches downsampled` == 0, + "Control", + ifelse( + imba_clus_merged$`Proportion downsampled` == 0, + "Ablated", + "Downsampled" + ) +) + +median_celltype_ari_results <- imba_clus_merged %>% + group_by(Method, type, `Downsampled celltypes`) %>% + summarize( + `Median celltype ARI` = median(`Celltype ARI Imbalanced`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +median_celltype_ari_results_vals_long <- reshape2::dcast( + median_celltype_ari_results, + formula = type + `Downsampled celltypes` ~ `Method`, + value.var = "Median celltype ARI" +) +median_celltype_ari_results_vals_long$type <- factor( + median_celltype_ari_results_vals_long$type, + levels = c("Control", "Downsampled", "Ablated") +) +median_celltype_ari_results_vals_long <- median_celltype_ari_results_vals_long[ + order( + median_celltype_ari_results_vals_long$type, + median_celltype_ari_results_vals_long$`Downsampled celltypes` + ), +] +rownames(median_celltype_ari_results_vals_long) <- c( + paste0("R_", seq(1, nrow(median_celltype_ari_results_vals_long))) +) +colnames(median_celltype_ari_results_vals_long)[1] <- c( + "Type" +) +median_celltype_ari_long_vals_only <- median_celltype_ari_results_vals_long[ + ,-c(1,2) +] + +median_celltype_ari_long_type <- median_celltype_ari_results_vals_long[ + ,1, drop = FALSE +] +median_celltype_ari_long_celltype <- median_celltype_ari_results_vals_long[ + ,2, drop = FALSE +] + +dark_2_cols = palette.colors(n = 8, "Dark2") +col_type = c( + "Control" = "forestgreen", + "Downsampled" = "darkorchid3", + "Ablated" = "firebrick2" +) +col_celltype = c( + "B cell" = dark_2_cols[1], + "CD14+ Monocyte" = dark_2_cols[2], + "CD4+ T cell" = dark_2_cols[3], + "CD8+ T cell" = dark_2_cols[4], + "FCGR3A+ Monocyte" = dark_2_cols[5], + "NK cell" = dark_2_cols[6], + "None" = "black" +) + +ht1 = Heatmap( + as.matrix(median_celltype_ari_long_vals_only), + name = "Median \ncell-type ARI", + width = unit(5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE +) +ht2 = Heatmap( + as.matrix(median_celltype_ari_long_type), + name = "Type", + col = col_type, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +ht3 = Heatmap( + as.matrix(median_celltype_ari_long_celltype), + name = "Affected cell-type", + col = col_celltype, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +celltype_ari_hm <- ht1 + ht2 + ht3 +CairoPDF( + "outs/control_w_cider/figures/17_celltype_ari_unscaled_ds_effects_heatmap_no_liger.pdf", + width = 8, + height = 6 +) +draw( + celltype_ari_hm, + column_title = "Integration method", + column_title_side = "bottom", + column_title_gp = gpar(fontsize = 14, fontface = "bold") +) +dev.off() + +# Repeat the same analysis above, but now with unscaled median batch ARI +median_batch_ari_results <- imba_clus_merged %>% + group_by(Method, type, `Downsampled celltypes`) %>% + summarize( + `Median batch ARI` = median(`Batch ARI`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +median_batch_ari_results_vals_long <- reshape2::dcast( + median_batch_ari_results, + formula = type + `Downsampled celltypes` ~ `Method`, + value.var = "Median batch ARI" +) +median_batch_ari_results_vals_long$type <- factor( + median_batch_ari_results_vals_long$type, + levels = c("Control", "Downsampled", "Ablated") +) +median_batch_ari_results_vals_long <- median_batch_ari_results_vals_long[ + order( + median_batch_ari_results_vals_long$type, + median_batch_ari_results_vals_long$`Downsampled celltypes` + ), +] +rownames(median_batch_ari_results_vals_long) <- c( + paste0("R_", seq(1, nrow(median_batch_ari_results_vals_long))) +) +colnames(median_batch_ari_results_vals_long)[1] <- c( + "Type" +) +median_batch_ari_long_vals_only <- median_batch_ari_results_vals_long[ + ,-c(1,2) +] +median_batch_ari_long_type <- median_batch_ari_results_vals_long[ + ,1, drop = FALSE +] +median_batch_ari_long_celltype <- median_batch_ari_results_vals_long[ + ,2, drop = FALSE +] + +dark_2_cols = palette.colors(n = 8, "Dark2") +col_type = c( + "Control" = "forestgreen", + "Downsampled" = "darkorchid3", + "Ablated" = "firebrick2" +) +col_celltype = c( + "B cell" = dark_2_cols[1], + "CD14+ Monocyte" = dark_2_cols[2], + "CD4+ T cell" = dark_2_cols[3], + "CD8+ T cell" = dark_2_cols[4], + "FCGR3A+ Monocyte" = dark_2_cols[5], + "NK cell" = dark_2_cols[6], + "None" = "black" +) + +ht1 = Heatmap( + as.matrix(median_batch_ari_long_vals_only), + name = "Median \nbatch ARI", + width = unit(5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE +) +ht2 = Heatmap( + as.matrix(median_batch_ari_long_type), + name = "Type", + col = col_type, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +ht3 = Heatmap( + as.matrix(median_batch_ari_long_celltype), + name = "Affected cell-type", + col = col_celltype, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +batch_ari_hm <- ht1 + ht2 + ht3 +CairoPDF( + "outs/control_w_cider/figures/17_batch_ari_unscaled_ds_effects_heatmap_no_liger.pdf", + width = 8, + height = 6 +) +draw( + batch_ari_hm, + column_title = "Integration method", + column_title_side = "bottom", + column_title_gp = gpar(fontsize = 14, fontface = "bold") +) +dev.off() + +# Repeat the same analysis above, but now with unscaled median balanced celltype ARI +median_balanced_ari_results <- imba_clus_merged %>% + group_by(Method, type, `Downsampled celltypes`) %>% + summarize( + `Median cell-type balanced ARI` = median(`Celltype ARI Balanced`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +median_balanced_ari_results_vals_long <- reshape2::dcast( + median_balanced_ari_results, + formula = type + `Downsampled celltypes` ~ `Method`, + value.var = "Median cell-type balanced ARI" +) +median_balanced_ari_results_vals_long$type <- factor( + median_balanced_ari_results_vals_long$type, + levels = c("Control", "Downsampled", "Ablated") +) +median_balanced_ari_results_vals_long <- median_balanced_ari_results_vals_long[ + order( + median_balanced_ari_results_vals_long$type, + median_balanced_ari_results_vals_long$`Downsampled celltypes` + ), +] +rownames(median_balanced_ari_results_vals_long) <- c( + paste0("R_", seq(1, nrow(median_balanced_ari_results_vals_long))) +) +colnames(median_balanced_ari_results_vals_long)[1] <- c( + "Type" +) +median_balanced_ari_long_vals_only <- median_balanced_ari_results_vals_long[ + ,-c(1,2) +] +median_balanced_ari_long_type <- median_balanced_ari_results_vals_long[ + ,1, drop = FALSE +] +median_balanced_ari_long_celltype <- median_balanced_ari_results_vals_long[ + ,2, drop = FALSE +] + +ht1 = Heatmap( + as.matrix(median_balanced_ari_long_vals_only), + name = "Median cell-type unscaled \nbalanced ARI", + width = unit(5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE +) +ht2 = Heatmap( + as.matrix(median_balanced_ari_long_type), + name = "Type", + col = col_type, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) +ht3 = Heatmap( + as.matrix(median_balanced_ari_long_celltype), + name = "Affected cell-type", + col = col_celltype, + width = unit(0.5, "cm"), + cluster_rows = FALSE, + cluster_columns = FALSE, + show_row_names = FALSE, + show_column_names = FALSE +) + +celltype_balanced_ari_hm <- ht1 + ht2 + ht3 +CairoPDF( + "outs/control_w_cider/figures/17_celltype_balanced_ari_unscaled_ds_effects_heatmap_no_liger.pdf", + width = 8, + height = 6 +) +draw( + celltype_balanced_ari_hm, + column_title = "Integration method", + column_title_side = "bottom", + column_title_gp = gpar(fontsize = 14, fontface = "bold") +) +dev.off() + +# Results of celltype downsampling and ablation on +# cluster-based classification accuracy + +# Merge imbalance and cluster classification results together +imba_clus_merged <- merge( + imba_concat, + clus_class_concat, + by = c( + "Number of batches downsampled", + "Number of celltypes downsampled", + "Proportion downsampled", + "Replicate" + ) +) +imba_clus_merged <- distinct(imba_clus_merged) + +# Subset for only cases where the celltype downsampled is equal to the +# celltype being classified +imba_clus_merged_celltype <- imba_clus_merged[ + imba_clus_merged$Celltype == imba_clus_merged$`Downsampled celltypes` | + imba_clus_merged$`Downsampled celltypes` %in% c("None") +] + +# Indicate which panels are control and which ones are ablations or downsampling +imba_clus_merged_celltype$type <- ifelse( + imba_clus_merged_celltype$`Number of batches downsampled` == 0, + "Control", + ifelse( + imba_clus_merged_celltype$`Proportion downsampled` == 0, + "Ablated", + "Downsampled" + ) +) + +# Format celltype names +imba_clus_merged_celltype$Celltype <- plyr::mapvalues( + imba_clus_merged_celltype$Celltype, + from = c( + "Monocyte_CD14", + "Monocyte_FCGR3A", + "CD4 T cell", + "CD8 T cell" + ), + to = c( + "CD14+ Monocyte", + "FCGR3A+ Monocyte", + "CD4+ T cell", + "CD8+ T cell" + ) +) + +# Create function to format facet labels (downsampled celltypes) +ds_celltype_labelled <- function(variable,value){ + return(paste0("Cell-type affected = ", value)) +} + +ggplot(data = imba_clus_merged_celltype, aes(x = `Method`, y = `F1-score`)) + + geom_boxplot( + aes( + fill = factor(`type`, levels = c("Control", "Downsampled", "Ablated")), + ), + notch = FALSE, + alpha = 0.8 + ) + + facet_wrap( + .~Celltype, + scales = "free_x", + labeller = ds_celltype_labelled, + ncol = 2 + ) + + labs( + fill = "Type", + x = "Method", + y = "Affected cell-type F1-classification score post-integration" + ) + + scale_fill_manual( + breaks = c("Control", "Downsampled", "Ablated"), + values = c("forestgreen", "darkorchid3", "firebrick2") + ) + + theme_few() + + theme(axis.title.x = element_text(size = 16)) + + theme(axis.title.y = element_text(size = 16)) + + theme(strip.text.x = element_text(size = 16)) + + theme(plot.title = element_text(size = 14)) + + theme(axis.text.x = element_text(size = 16)) + + theme(axis.text.y = element_text(size = 16)) + + theme(legend.title = element_text(size = 16)) + + theme(legend.text = element_text(size = 16)) +ggsave( + "outs/control_w_cider/figures/17_pbmc_ds_ablate_allmethod_clus_f1_score_no_liger.pdf", + width = 16, + height = 14, + device = cairo_pdf +) + +ggplot(data = imba_clus_merged_celltype, aes(x = `Method`, y = `F1-score`)) + + geom_boxplot( + aes( + fill = factor(`type`, levels = c("Control", "Downsampled", "Ablated")), + ), + notch = FALSE, + alpha = 0.8 + ) + + ylim(0, 1) + + facet_wrap( + .~Celltype, + scales = "free_x", + labeller = ds_celltype_labelled, + ncol = 2 + ) + + labs( + fill = "Type", + x = "Method", + y = "Affected cell-type F1-classification score post-integration" + ) + + scale_fill_manual( + breaks = c("Control", "Downsampled", "Ablated"), + values = c("forestgreen", "darkorchid3", "firebrick2") + ) + + theme_few() + + theme(axis.title.x = element_text(size = 16)) + + theme(axis.title.y = element_text(size = 16)) + + theme(strip.text.x = element_text(size = 16)) + + theme(plot.title = element_text(size = 14)) + + theme(axis.text.x = element_text(size = 16)) + + theme(axis.text.y = element_text(size = 16)) + + theme(legend.title = element_text(size = 16)) + + theme(legend.text = element_text(size = 16)) +ggsave( + "outs/control_w_cider/figures/17_pbmc_ds_ablate_allmethod_clus_f1_score_no_liger_0_1_y.pdf", + width = 16, + height = 14, + device = cairo_pdf +) + +# Get and save the standard deviation of the median values of F1, based on +# method utilized, celltype downsampled, and subset utilized +imba_clus_merged_celltype_medians <- imba_clus_merged_celltype %>% + group_by(Method, Celltype, type) %>% + summarize( + `Median F1-score per subset` = median(`F1-score`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +imba_clus_merged_celltype_medians_stdev <- imba_clus_merged_celltype_medians %>% + group_by(Celltype) %>% + summarize( + `Stdev Method Median F1-score per subset` = sd( + `Median F1-score per subset`, + na.rm = FALSE + ), + .groups = "keep" + ) %>% + as.data.frame() + +colnames(imba_clus_merged_celltype_medians_stdev) <- c( + "Cell-type", + "Standard deviation of medians for F1-score, across methods, replicates, and experiment types" +) +fwrite( + imba_clus_merged_celltype_medians_stdev, + "outs/control_w_cider/results/17_baseline_pbmc_clus_class_f1_score_stdevs_per_celltype.tsv", + sep = "\t", + row.names = FALSE, + col.names = TRUE +) + +# Results of celltype downsampling and ablation on +# KNN-cluster-based classification accuracy +# The difference between this analysis is that the previous used logistic +# regression on the cluster labels while this one uses KNN on the cluster labels + +# Merge imbalance and knn cluster classification results together +imba_knn_clus_merged <- merge( + imba_concat, + knn_clus_class_concat, + by = c( + "Number of batches downsampled", + "Number of celltypes downsampled", + "Proportion downsampled", + "Replicate" + ) +) +imba_knn_clus_merged <- distinct(imba_knn_clus_merged) + +# Subset for only cases where the celltype downsampled is equal to the +# celltype being classified +imba_knn_clus_merged_celltype <- imba_knn_clus_merged[ + imba_knn_clus_merged$Celltype == imba_knn_clus_merged$`Downsampled celltypes` | + imba_knn_clus_merged$`Downsampled celltypes` %in% c("None") +] + +# Indicate which panels are control and which ones are ablations or downsampling +imba_knn_clus_merged_celltype$type <- ifelse( + imba_knn_clus_merged_celltype$`Number of batches downsampled` == 0, + "Control", + ifelse( + imba_knn_clus_merged_celltype$`Proportion downsampled` == 0, + "Ablated", + "Downsampled" + ) +) + +# Format celltype names +imba_knn_clus_merged_celltype$Celltype <- plyr::mapvalues( + imba_knn_clus_merged_celltype$Celltype, + from = c( + "Monocyte_CD14", + "Monocyte_FCGR3A", + "CD4 T cell", + "CD8 T cell" + ), + to = c( + "CD14+ Monocyte", + "FCGR3A+ Monocyte", + "CD4+ T cell", + "CD8+ T cell" + ) +) + +# Create function to format facet labels (downsampled celltypes) +ds_celltype_labelled <- function(variable,value){ + return(paste0("Cell-type affected = ", value)) +} + +ggplot(data = imba_knn_clus_merged_celltype, aes(x = `Method`, y = `F1-score`)) + + geom_boxplot( + aes( + fill = factor(`type`, levels = c("Control", "Downsampled", "Ablated")), + ), + notch = FALSE, + alpha = 0.8 + ) + + facet_wrap( + .~Celltype, + scales = "free_x", + labeller = ds_celltype_labelled, + ncol = 2 + ) + + labs( + fill = "Type", + x = "Method", + y = "Affected cell-type F1-classification score post-integration" + ) + + scale_fill_manual( + breaks = c("Control", "Downsampled", "Ablated"), + values = c("forestgreen", "darkorchid3", "firebrick2") + ) + + theme_few() + + theme(axis.title.x = element_text(size = 16)) + + theme(axis.title.y = element_text(size = 16)) + + theme(strip.text.x = element_text(size = 16)) + + theme(plot.title = element_text(size = 14)) + + theme(axis.text.x = element_text(size = 16)) + + theme(axis.text.y = element_text(size = 16)) + + theme(legend.title = element_text(size = 16)) + + theme(legend.text = element_text(size = 16)) +ggsave( + "outs/control_w_cider/figures/17_pbmc_ds_ablate_allmethod_knn_clus_f1_score_no_liger.pdf", + width = 16, + height = 14, + device = cairo_pdf +) + +ggplot(data = imba_knn_clus_merged_celltype, aes(x = `Method`, y = `F1-score`)) + + geom_boxplot( + aes( + fill = factor(`type`, levels = c("Control", "Downsampled", "Ablated")), + ), + notch = FALSE, + alpha = 0.8 + ) + + ylim(0, 1) + + facet_wrap( + .~Celltype, + scales = "free_x", + labeller = ds_celltype_labelled, + ncol = 2 + ) + + labs( + fill = "Type", + x = "Method", + y = "Affected cell-type F1-classification score post-integration" + ) + + scale_fill_manual( + breaks = c("Control", "Downsampled", "Ablated"), + values = c("forestgreen", "darkorchid3", "firebrick2") + ) + + theme_few() + + theme(axis.title.x = element_text(size = 16)) + + theme(axis.title.y = element_text(size = 16)) + + theme(strip.text.x = element_text(size = 16)) + + theme(plot.title = element_text(size = 14)) + + theme(axis.text.x = element_text(size = 16)) + + theme(axis.text.y = element_text(size = 16)) + + theme(legend.title = element_text(size = 16)) + + theme(legend.text = element_text(size = 16)) +ggsave( + "outs/control_w_cider/figures/17_pbmc_ds_ablate_allmethod_knn_clus_f1_score_no_liger_0_1_y.pdf", + width = 16, + height = 14, + device = cairo_pdf +) + +# Get and save the standard deviation of the median values of F1, based on +# method utilized, celltype downsampled, and subset utilized +imba_knn_clus_merged_celltype_medians <- imba_knn_clus_merged_celltype %>% + group_by(Method, Celltype, type) %>% + summarize( + `Median F1-score per subset` = median(`F1-score`, na.rm = FALSE), + .groups = "keep" + ) %>% + as.data.frame() + +imba_knn_clus_merged_celltype_medians_stdev <- imba_knn_clus_merged_celltype_medians %>% + group_by(Celltype) %>% + summarize( + `Stdev Method Median F1-score per subset` = sd( + `Median F1-score per subset`, + na.rm = FALSE + ), + .groups = "keep" + ) %>% + as.data.frame() + +colnames(imba_knn_clus_merged_celltype_medians_stdev) <- c( + "Cell-type", + "Standard deviation of medians for F1-score, across methods, replicates, and experiment types" +) +fwrite( + imba_knn_clus_merged_celltype_medians_stdev, + "outs/control_w_cider/results/17_baseline_pbmc_knn_clus_f1_score_stdevs_per_celltype.tsv", + sep = "\t", + row.names = FALSE, + col.names = TRUE +) \ No newline at end of file diff --git a/workflow/analysis/R/21_PBMC_perturbation_umap_plots.R b/workflow/analysis/R/21_PBMC_perturbation_umap_plots.R new file mode 100644 index 0000000..15549dd --- /dev/null +++ b/workflow/analysis/R/21_PBMC_perturbation_umap_plots.R @@ -0,0 +1,206 @@ +library(data.table) +library(tidyverse) +library(reshape2) +library(ggplot2) +library(ggthemes) +library(ggExtra) +library(ggpubr) +library(dotwhisker) +library(Seurat) +library(SeuratDisk) +library(ComplexHeatmap) +library(circlize) +library(RColorBrewer) +library(Cairo) + +# Helper functions +`%ni%` <- Negate(`%in%`) + +# Change to results dir for uamp results data +setwd("../../../results/umap/") + +# Load color palette +kev_palette <- c("dodgerblue2", + "#E31A1C", + "green4", + "#6A3D9A", + "#FF7F00", + "black", + "gold1", + "skyblue2", + "#FB9A99", + "palegreen2", + "#CAB2D6", + "#FDBF6F", + "gray70", + "khaki2", + "maroon", + "orchid1", + "deeppink1", + "blue1", + "steelblue4", + "darkturquoise", + "green1", + "yellow4", + "yellow3", + "darkorange4", + "brown") + +##### Analysis of PBMC 2 batch balanced data - baseline ##### + +# Load in the umap plot results +setwd("umap_plots/") +umap_files <- list.files() +umap_files <- grep( + ".tsv", + umap_files, + value = TRUE +) +umap_files <- grep( + "pbmc_2_batch_base_balanced", + umap_files, + value = TRUE +) +umap_loaded <- lapply(umap_files, fread) +umap_names <- str_split_fixed(umap_files, fixed(".tsv"), 2)[,1] +names(umap_loaded) <- umap_names + +setwd("../../..") + +# Create directory for umap results if it doesn't exist +if (!dir.exists("outs/umap/results")) { + dir.create("outs/umap/results", recursive = TRUE) +} +if (!dir.exists("outs/umap/figures")) { + dir.create("outs/umap/figures") +} + +# Create function to loop over the umap files and return the results +umap_plot <- function(df, save_prefix) { + # Format celltype names + df$Clustering <- plyr::mapvalues( + df$Clustering, + from = c( + "Monocyte_CD14", + "Monocyte_FCGR3A", + "CD4 T cell", + "CD8 T cell" + ), + to = c( + "CD14+ Monocyte", + "FCGR3A+ Monocyte", + "CD4+ T cell", + "CD8+ T cell" + ) + ) + + # Format batch names + df$Clustering <- plyr::mapvalues( + df$Clustering, + from = c( + "batch_1", + "batch_2" + ), + to = c( + "Batch 1", + "Batch 2" + ) + ) + + unique_cluster_len <- length(unique(df$Clustering)) + if (unique_cluster_len > 8) { + ggplot(data = df, aes(x = `UMAP 1`, y = `UMAP 2`)) + + geom_point( + aes( + color = factor( + as.numeric(Clustering), + levels = sort(as.numeric(unique(df$Clustering))) + ) + ), + size = 0.25 + ) + + facet_wrap( + .~Subset, + scales = "free" + ) + + labs( + color = "", + x = "UMAP 1", + y = "UMAP 2" + ) + + scale_color_manual( + name = "", + values = kev_palette[1:unique_cluster_len] + ) + + guides(color = guide_legend(override.aes = list(size=2))) + + theme_few() + + theme(axis.title.x = element_text(size = 16)) + + theme(axis.title.y = element_text(size = 16)) + + theme(strip.text.x = element_text(size = 16)) + + theme(plot.title = element_text(size = 14)) + + theme(axis.text.x = element_text(size = 16)) + + theme(axis.text.y = element_text(size = 16)) + + theme(legend.title = element_text(size = 16)) + + theme(legend.text = element_text(size = 16)) + ggsave( + paste0( + "outs/umap/figures/", + save_prefix, + ".pdf" + ), + width = 16, + height = 8, + device = cairo_pdf + ) + } else { + if (any(grepl("Batch", df$Clustering))) { + pal = "Set1" + } else { + pal = "Dark2" + } + ggplot(data = df, aes(x = `UMAP 1`, y = `UMAP 2`)) + + geom_point( + aes( + color = factor(Clustering), + ), + size = 0.5 + ) + + facet_wrap( + .~Subset, + scales = "free" + ) + + labs( + color = "", + x = "UMAP 1", + y = "UMAP 2" + ) + + guides(color = guide_legend(override.aes = list(size=2))) + + scale_color_brewer(palette = pal) + + theme_few() + + theme(axis.title.x = element_text(size = 16)) + + theme(axis.title.y = element_text(size = 16)) + + theme(strip.text.x = element_text(size = 16)) + + theme(plot.title = element_text(size = 14)) + + theme(axis.text.x = element_text(size = 16)) + + theme(axis.text.y = element_text(size = 16)) + + theme(legend.title = element_text(size = 16)) + + theme(legend.text = element_text(size = 16)) + ggsave( + paste0( + "outs/umap/figures/", + save_prefix, + ".pdf" + ), + width = 16, + height = 8, + device = cairo_pdf + ) + } +} + +# Iterate over the umap files and names and save the results +mapply( + umap_plot, + df = umap_loaded, + save_prefix = umap_names +) diff --git a/workflow/configs/config_control_cider.json b/workflow/configs/config_control_cider.json new file mode 100644 index 0000000..177e292 --- /dev/null +++ b/workflow/configs/config_control_cider.json @@ -0,0 +1,22 @@ +{ + "config_name": "control_w_cider", + "int_datasets": { + "pbmc_2_batch_base_balanced" : { + "data_folder": "pbmc_2_batch_base_balanced", + "ds_celltypes": [1], + "ds_proportions": [0.1, 0], + "num_batches": [0, 1], + "repetitions": 200 + }, + "pbmc_2_batch_hierarchical_balanced": { + "data_folder": "pbmc_2_batch_hierarchical_balanced", + "ds_celltypes": [1], + "ds_proportions": [0.1, 0], + "num_batches": [0, 1], + "repetitions": 200 + } + }, + "int_ti_datasets": {}, + "query_to_reference": "No", + "celltype_list": "No" +} \ No newline at end of file diff --git a/workflow/configs/config_umap.json b/workflow/configs/config_umap.json new file mode 100644 index 0000000..b3a8d9a --- /dev/null +++ b/workflow/configs/config_umap.json @@ -0,0 +1,50 @@ +{ + "config_name": "umap", + "int_datasets": { + "pbmc_2_batch_base_balanced" : { + "data_folder": "pbmc_2_batch_base_balanced", + "ds_celltypes": [1], + "ds_proportions": [0.1, 0], + "num_batches": [0, 1], + "repetitions": 1 + }, + "pbmc_2_batch_hierarchical_balanced": { + "data_folder": "pbmc_2_batch_hierarchical_balanced", + "ds_celltypes": [1], + "ds_proportions": [0.1, 0], + "num_batches": [0, 1], + "repetitions": 1 + }, + "pbmc_2_batch" : { + "data_folder": "pbmc_2_batch", + "ds_celltypes": [0], + "ds_proportions": [0], + "num_batches": [0], + "repetitions": 1 + }, + "pbmc_4_batch" : { + "data_folder": "pbmc_4_batch", + "ds_celltypes": [0], + "ds_proportions": [0], + "num_batches": [0], + "repetitions": 1 + }, + "mouse_hindbrain_6_batch": { + "data_folder": "mouse_hindbrain_6_batch", + "ds_celltypes": [0], + "ds_proportions": [0], + "num_batches": [0], + "repetitions": 1 + }, + "peng_pdac_8_batch": { + "data_folder": "peng_pdac_tumor_annot_8_batch", + "ds_celltypes": [0], + "ds_proportions": [0], + "num_batches": [0], + "repetitions": 1 + } + }, + "int_ti_datasets": {}, + "query_to_reference": "No", + "celltype_list": "No" +} \ No newline at end of file diff --git a/workflow/scripts/R/install_packages.R b/workflow/scripts/R/install_packages.R new file mode 100644 index 0000000..4be001b --- /dev/null +++ b/workflow/scripts/R/install_packages.R @@ -0,0 +1,19 @@ +# Install the required packages +install.packages("BiocManager", repos = "http://cran.us.r-project.org") +library(BiocManager) +BiocManager::install("limma") +install.packages( + "https://cran.r-project.org/src/contrib/Archive/locfit/locfit_1.5-9.4.tar.gz", + repos = NULL, + type = "source" +) +BiocManager::install("edgeR") +install.packages("CIDER", repos = "http://cran.us.r-project.org") + +# Stop if CIDER is not installed successfully +if (!"CIDER" %in% installed.packages()[, "Package"]) { + stop("Package CIDER not installed successfully.") +} + +# Output a dummy file to indicate that the script has completed +file.create("../results/install_confirmation.txt") diff --git a/workflow/scripts/python/cluster_classification.py b/workflow/scripts/python/cluster_classification.py new file mode 100644 index 0000000..8c2d974 --- /dev/null +++ b/workflow/scripts/python/cluster_classification.py @@ -0,0 +1,135 @@ +import argparse + +import scanpy as sc +import anndata as ann +import numpy as np +import pandas as pd +from sklearn.model_selection import train_test_split +from sklearn.metrics import classification_report +from sklearn.linear_model import LogisticRegression + +def main(h5ad_loc, save_loc, dataset_name, rep): + # Load h5ad file + adata = sc.read_h5ad(h5ad_loc) + + # Extract summary statistics from h5ad file + num_batches_ds = adata.uns["downsampling_stats"]["num_batches"] + num_celltypes_ds = adata.uns["downsampling_stats"]["num_celltypes_downsampled"] + prop_ds = adata.uns["downsampling_stats"]["proportion_downsampled"] + + # Subset h5ad based on batch-correction method used + adata_method_sub = [] + methods = ["harmony", "scvi", "scanorama", "bbknn", "seurat", "liger", "cider"] + for method in methods: + adata_sub = adata[adata.obs["integration_method"] == method] + adata_method_sub.append( + adata_sub + ) + + # Determine cluster-based accuracy for each batch-correction method + precision_scores = [] + recall_scores = [] + f1_scores = [] + supports = [] + celltypes = [] + for adata_sub in adata_method_sub: + # Split testing and training data in stratified manner (70/30) + X = adata_sub.obs["leiden"].__array__() + y = adata_sub.obs["celltype"].__array__() + X_train, X_test, y_train, y_test = train_test_split( + X, y, stratify=y, test_size=0.7, random_state=42 + ) + + # Train a logistic regression classifier to predict celltype based on cluster + X_train_num = np.array( + [int(x) for x in X_train] + ) + X_test_num = np.array( + [int(x) for x in X_test] + ) + # Reshape X test and train to be 2D arrays + X_train_num = X_train_num.reshape(-1, 1) + X_test_num = X_test_num.reshape(-1, 1) + + # Train logistic regression classifier and predict celltypes + clf = LogisticRegression( + random_state=42, max_iter=1000 + ).fit(X_train_num, y_train) + y_pred = clf.predict(X_test_num) + + # Get classification report and subset for only relevant columns + class_report_dict = classification_report( + y_test, y_pred, output_dict=True + ) + class_report_df = pd.DataFrame(class_report_dict) + class_report_df = class_report_df.iloc[:, :-3] + + # Append appropriate values to lists + precision_scores.append(class_report_df.loc["precision"].values) + recall_scores.append(class_report_df.loc["recall"].values) + f1_scores.append(class_report_df.loc["f1-score"].values) + supports.append(class_report_df.loc["support"].values) + celltypes.append(class_report_df.columns.values) + + # Repeat method values to have same length as scores (one for each celltype) + methods_repeat = np.repeat(methods, len(precision_scores[0])) + + # Concatenate scores and celltypes + precision_scores_concat = np.concatenate(precision_scores) + recall_scores_concat = np.concatenate(recall_scores) + f1_scores_concat = np.concatenate(f1_scores) + supports_concat = np.concatenate(supports) + celltypes_concat = np.concatenate(celltypes) + + # Create summary dataframe for classification statistics and save + classification_summary_df = pd.DataFrame({ + "Dataset": dataset_name, + "Number of batches downsampled": num_batches_ds, + "Number of celltypes downsampled": num_celltypes_ds, + "Proportion downsampled": prop_ds, + "Replicate": rep, + "Method": methods_repeat, + "Celltype": celltypes_concat, + "Precision": precision_scores_concat, + "Recall": recall_scores_concat, + "F1-score": f1_scores_concat, + "Support": supports_concat, + "Mean LogReg F1-score": np.mean(f1_scores_concat) + }) + classification_summary_df.to_csv( + save_loc, + index=False, + sep="\t" + ) + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description = "Input and output files for cluster classification summary" + ) + parser.add_argument( + "--infile", + type = str, + help = "Path of integrated h5ad file" + ) + parser.add_argument( + "--outfile", + type = str, + help = "Filepath for saving cluster classification results of integrated h5ad file" + ) + parser.add_argument( + "--dataset", + type = str, + help = "Name of dataset" + ) + parser.add_argument( + "--rep", + type = int, + help = "Repetition number" + ) + args = parser.parse_args() + main( + h5ad_loc = args.infile, + save_loc = args.outfile, + dataset_name = args.dataset, + rep = args.rep + ) \ No newline at end of file diff --git a/workflow/scripts/python/cluster_knn_classification.py b/workflow/scripts/python/cluster_knn_classification.py new file mode 100644 index 0000000..ef9e8b9 --- /dev/null +++ b/workflow/scripts/python/cluster_knn_classification.py @@ -0,0 +1,135 @@ +import argparse + +import scanpy as sc +import anndata as ann +import numpy as np +import pandas as pd +from sklearn.model_selection import train_test_split +from sklearn.metrics import classification_report +from sklearn.neighbors import KNeighborsClassifier +from sklearn.preprocessing import OneHotEncoder + +def main(h5ad_loc, save_loc, dataset_name, rep): + # Load h5ad file + adata = sc.read_h5ad(h5ad_loc) + + # Extract summary statistics from h5ad file + num_batches_ds = adata.uns["downsampling_stats"]["num_batches"] + num_celltypes_ds = adata.uns["downsampling_stats"]["num_celltypes_downsampled"] + prop_ds = adata.uns["downsampling_stats"]["proportion_downsampled"] + + # Subset h5ad based on batch-correction method used + adata_method_sub = [] + methods = ["harmony", "scvi", "scanorama", "bbknn", "seurat", "liger", "cider"] + for method in methods: + adata_sub = adata[adata.obs["integration_method"] == method] + adata_method_sub.append( + adata_sub + ) + + # Determine KNN-cluster-based accuracy for each batch-correction method + precision_scores = [] + recall_scores = [] + f1_scores = [] + supports = [] + celltypes = [] + for adata_sub in adata_method_sub: + # Split testing and training data in stratified manner (70/30) + X = adata_sub.obs["leiden"].__array__() + y = adata_sub.obs["celltype"].__array__() + X_train, X_test, y_train, y_test = train_test_split( + X, y, stratify=y, test_size=0.7, random_state=42 + ) + + # One-hot encode X train and test + enc = OneHotEncoder(handle_unknown="ignore") # ignore unknown categories in test set + X_train_reshaped = X_train.reshape(-1, 1) + X_test_reshaped = X_test.reshape(-1, 1) + enc.fit(X_train_reshaped) + X_train = enc.transform(X_train_reshaped).toarray() + X_test = enc.transform(X_test_reshaped).toarray() + + # Train k-nearest neighbors classifier with k=15 and predict on test data + knn = KNeighborsClassifier( + n_neighbors=15 + ) + knn.fit(X_train, y_train) + y_pred = knn.predict(X_test) + + # Get classification report and subset for only relevant columns + class_report_dict = classification_report( + y_test, y_pred, output_dict=True + ) + class_report_df = pd.DataFrame(class_report_dict) + class_report_df = class_report_df.iloc[:, :-3] + + # Append appropriate values to lists + precision_scores.append(class_report_df.loc["precision"].values) + recall_scores.append(class_report_df.loc["recall"].values) + f1_scores.append(class_report_df.loc["f1-score"].values) + supports.append(class_report_df.loc["support"].values) + celltypes.append(class_report_df.columns.values) + + # Repeat method values to have same length as scores (one for each celltype) + methods_repeat = np.repeat(methods, len(precision_scores[0])) + + # Concatenate scores and celltypes + precision_scores_concat = np.concatenate(precision_scores) + recall_scores_concat = np.concatenate(recall_scores) + f1_scores_concat = np.concatenate(f1_scores) + supports_concat = np.concatenate(supports) + celltypes_concat = np.concatenate(celltypes) + + # Create summary dataframe for classification statistics and save + classification_summary_df = pd.DataFrame({ + "Dataset": dataset_name, + "Number of batches downsampled": num_batches_ds, + "Number of celltypes downsampled": num_celltypes_ds, + "Proportion downsampled": prop_ds, + "Replicate": rep, + "Method": methods_repeat, + "Celltype": celltypes_concat, + "Precision": precision_scores_concat, + "Recall": recall_scores_concat, + "F1-score": f1_scores_concat, + "Support": supports_concat, + "Mean KNN F1-score": np.mean(f1_scores_concat) + }) + classification_summary_df.to_csv( + save_loc, + index=False, + sep="\t" + ) + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description = "Input and output files for knn-cluster classification summary" + ) + parser.add_argument( + "--infile", + type = str, + help = "Path of integrated h5ad file" + ) + parser.add_argument( + "--outfile", + type = str, + help = "Filepath for saving cluster classification results of integrated h5ad file" + ) + parser.add_argument( + "--dataset", + type = str, + help = "Name of dataset" + ) + parser.add_argument( + "--rep", + type = int, + help = "Repetition number" + ) + args = parser.parse_args() + main( + h5ad_loc = args.infile, + save_loc = args.outfile, + dataset_name = args.dataset, + rep = args.rep + ) + \ No newline at end of file diff --git a/workflow/scripts/python/umap_plots.py b/workflow/scripts/python/umap_plots.py new file mode 100644 index 0000000..d5b7eaf --- /dev/null +++ b/workflow/scripts/python/umap_plots.py @@ -0,0 +1,129 @@ +import argparse +import os +import sys +sys.path.append("src/python/") + +import numpy as np +import pandas as pd +import anndata as ann +import scanpy as sc + +from utils import Umap + +def main(h5ad_loc, leiden_save_loc, celltype_save_loc, batch_save_loc): + # Load h5ad file + adata = sc.read_h5ad(h5ad_loc) + + # Get the umap coordinates for all the methods and create dictionary object + methods = ["bbknn", "harmony", "scanorama", "scvi", "seurat"] + umap_dict = {} + for method in methods: + adata_sub = adata[adata.obs["integration_method"] == method] + umap_dict[method] = adata_sub.obsm["X_umap"].__array__() + + # Get the leiden clustering for all the methods and create dictionary object + leiden_dict = {} + for method in methods: + adata_sub = adata[adata.obs["integration_method"] == method] + leiden_dict[method] = adata_sub.obs["leiden"].__array__() + + # Create a cell type dictionary object + celltype_dict = {} + for method in methods: + adata_sub = adata[adata.obs["integration_method"] == method] + celltype_dict[method] = adata_sub.obs["celltype"].__array__() + + # Create a batch dictionary object + batch_dict = {} + for method in methods: + adata_sub = adata[adata.obs["integration_method"] == method] + batch_dict[method] = adata_sub.obs["batch"].__array__() + + # Create a umap object for each subset of information + umap_leiden = Umap( + coords = umap_dict, + clustering = leiden_dict, + subset_name = "Clustering", + ) + umap_celltype = Umap( + coords = umap_dict, + clustering = celltype_dict, + subset_name = "Cell-type", + ) + umap_batch = Umap( + coords = umap_dict, + clustering = batch_dict, + subset_name = "Batch", + ) + + # Plot each of the umap objects + umap_leiden.umap_df() + umap_leiden.umap_plot(show_plot=True) + umap_leiden.save_umap( + save_dir=leiden_save_loc, + dpi=300 + ) + + umap_celltype.umap_df() + umap_celltype.umap_plot(show_plot=True) + umap_celltype.save_umap( + save_dir=celltype_save_loc, + dpi=300 + ) + + umap_batch.umap_df() + umap_batch.umap_plot(show_plot=True) + umap_batch.save_umap( + save_dir=batch_save_loc, + dpi=300 + ) + + # Save the umap dataframes as auxiliary files + umap_leiden_df = umap_leiden.umap_concat + umap_celltype_df = umap_celltype.umap_concat + umap_batch_df = umap_batch.umap_concat + + umap_leiden_df.to_csv( + os.path.splitext(leiden_save_loc)[0] + ".tsv", + sep = "\t" + ) + umap_celltype_df.to_csv( + os.path.splitext(celltype_save_loc)[0] + ".tsv", + sep = "\t" + ) + umap_batch_df.to_csv( + os.path.splitext(batch_save_loc)[0] + ".tsv", + sep = "\t" + ) + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + description = "Input and output files for UMAP plot generation" + ) + parser.add_argument( + "--infile", + type = str, + help = "Path of integrated h5ad file" + ) + parser.add_argument( + "--leiden_plot_loc", + type = str, + help = "Filepath for saving leiden overlayed UMAP results" + ) + parser.add_argument( + "--celltype_plot_loc", + type = str, + help = "Filepath for saving celltype overlayed UMAP results" + ) + parser.add_argument( + "--batch_plot_loc", + type = str, + help = "Filepath for saving batch overlayed UMAP results" + ) + args = parser.parse_args() + main( + h5ad_loc = args.infile, + leiden_save_loc = args.leiden_plot_loc, + celltype_save_loc = args.celltype_plot_loc, + batch_save_loc = args.batch_plot_loc + ) \ No newline at end of file diff --git a/workflow/src/R/cider_integrate.R b/workflow/src/R/cider_integrate.R new file mode 100644 index 0000000..0d78420 --- /dev/null +++ b/workflow/src/R/cider_integrate.R @@ -0,0 +1,53 @@ +library(CIDER) +library(Seurat) +library(SeuratDisk) +library(reticulate) +library(parallel) +library(doParallel) +library(plyr) + +# Read in matrix for full data, including last batch column +args <- commandArgs(trailingOnly = TRUE) +file <- args[1] +filename <- args[2] + +# Load anndata and scanpy +ad <- import("anndata") +sc <- import("scanpy") + +# Load h5ad object through reticulate and create Seurat object +temp_adata <- ad$read_h5ad(file) +exprs <- t(temp_adata$X$todense()) +colnames(exprs) <- temp_adata$obs_names$to_list() +rownames(exprs) <- temp_adata$var_names$to_list() +seur_obj <- CreateSeuratObject(exprs) +seur_obj <- SetAssayData(seur_obj, "data", exprs) +seur_obj <- AddMetaData(seur_obj, temp_adata$obs) + +# Format the batch names +batch_names <- seur_obj$batch +batch_names_unique <- unique(batch_names) +batch_names_encoded <- as.numeric(factor(batch_names)) +seur_obj$Batch <- batch_names_encoded + +# Use the high-level CIDER pipeline +seur_obj <- initialClustering(seur_obj, nfeatures = 2500, resolution = 1.0, dims = 1:20) +ider <- getIDEr(seur_obj, downsampling.size = 35, use.parallel = TRUE, n.cores = 8, verbose = FALSE) +seur_obj <- finalClustering(seur_obj, ider, cutree.h = 0.35) # final clustering + +# Extract the cluster labels and append to a dataframe +cluster_labels <- as.character(seur_obj$CIDER_cluster) +batch_labels <- as.character(seur_obj$batch) +cluster_df <- data.frame( + "batch" = batch_labels, + "cider_cluster" = cluster_labels +) + +# Write out cluster labels to file +write.table( + cluster_df, + paste0("./tmp/", filename, "_cider_out.txt"), + sep = "\t", + row.names = FALSE, + col.names = TRUE +) diff --git a/workflow/src/python/utils/cider_integrate.py b/workflow/src/python/utils/cider_integrate.py new file mode 100644 index 0000000..13e9c11 --- /dev/null +++ b/workflow/src/python/utils/cider_integrate.py @@ -0,0 +1,104 @@ +import os +import sys +import subprocess +import uuid + +from natsort import natsorted +import numpy as np +import pandas as pd +import scipy as sp +import anndata as ann +import scanpy as sc + +class CiderIntegrate: + """ + Class for interpolating between the Integration class and R-script + based integration of RNA-seq batches using the CIDER package. + As CIDER integration returns a set of cluster annotations, the + anndata object has the cluster annotations appended and returned + """ + def __init__(self, adata): + """ + Args: + adata (object): An instance of an anndata class corresponding to the CIDER + subset from the Integration class. + """ + self.adata = adata.copy() + + def _format(self): + # Append a column on gene names + self.adata.var["gene"] = self.adata.var_names + # Remove layers and raw from AnnData object (avoid conflicts with h5seurat) + self.adata.layers = None + self.adata.raw = None + # Add 'Batch' column to obs as CIDER requires this + self.adata.obs["Batch"] = self.adata.obs["batch"].values + + def _output_temp_h5ad(self): + # Check if temp exists, if not, make dir + if not os.path.exists("tmp"): + os.makedirs("tmp") + + # Output temporary file with data + self.filename = ''.join(str(uuid.uuid4()).split("-")) + self.file = "{filename}.h5ad".format(filename = self.filename) + self.adata.write_h5ad(os.path.join("tmp", self.file)) + + def _cider_integrate(self): + # Call subprocess and call R script + tempfile_script = \ + "Rscript src/R/cider_integrate.R tmp/{tempfile} {tempfile_name} --verbose".format( + tempfile = self.file, + tempfile_name = self.filename + ) + + self.sp_integrate = subprocess.run(tempfile_script, shell = True, text = True, capture_output = True) + if self.sp_integrate.returncode != 0: + raise Exception( + "Subprocess call returned nonzero exit code - call: {call} \n Output: {output}".format( + call = self.sp_integrate.stderr, + output = self.sp_integrate.stdout + ) + ) + + def _return_integrated(self): + # Get cider output file + self.cider_outfile = "{filename}_cider_out.txt".format(filename = self.filename) + + # Read in as pandas dataframe + integrated_cluster_annotations = pd.read_csv(os.path.join("tmp", self.cider_outfile), sep = "\t") + + # Append the cluster annotations to the original AnnData object - store as leiden + cider_annotations = integrated_cluster_annotations["cider_cluster"].values + self.adata.obs["leiden"] = pd.Categorical( + values=cider_annotations.astype('U'), + categories=natsorted(map(str, np.unique(cider_annotations))), + ) + + # Return integrated AnnData object + return self.adata + + def _clean_files(self): + # Remove temporary python and cider files + tmp_files = os.listdir("tmp") + tmp_files_instance = [f for f in tmp_files if self.filename in f] + for f in tmp_files_instance: + os.remove(os.path.join("tmp", f)) + + # Check if all files related to the filename are removed from folder + tmp_files = os.listdir("tmp") + tmp_files_instance = [f for f in tmp_files if self.filename in f] + if len(tmp_files_instance) > 0: + raise Exception( + "Temporary file cleanup incomplete - files remain in folder" + ) + + def integrate(self): + # Perform workflow and return integrated anndata object + self._format() + self._output_temp_h5ad() + self._cider_integrate() + integrated_adata = self._return_integrated() + self._clean_files() + + return integrated_adata \ No newline at end of file diff --git a/workflow/src/python/utils/umap.py b/workflow/src/python/utils/umap.py new file mode 100644 index 0000000..a7ea270 --- /dev/null +++ b/workflow/src/python/utils/umap.py @@ -0,0 +1,122 @@ +import pandas as pd +import numpy as np +import seaborn as sns +import colorcet as cc +from natsort import natsorted + + +class Umap: + """Class for plotting results of integration experiments""" + + def __init__(self, coords, clustering, subset_name = None): + """ + Args: + coords (dictionary): coordinates of umap in numpy format where + keys correspond to following integration methods - + 'bbknn', 'harmony', 'scanorama', 'seurat', and 'scvi'. + clustering (dictionary): leiden or celltype clustering in numpy + format of integrated where keys correspond to following + integration methods - 'bbknn', 'harmony', 'scanorama', + 'seurat', and 'scvi'. + subset_name (string): name of subset being utilized for clustering + comparisons (e.g. batch, celltype). + """ + self.clustering_harmony = clustering.get("harmony") + self.clustering_scvi = clustering.get("scvi") + self.clustering_bbknn = clustering.get("bbknn") + self.clustering_scanorama = clustering.get("scanorama") + self.clustering_seurat = clustering.get("seurat") + + self.umap_harmony = coords.get("harmony") + self.umap_scvi = coords.get("scvi") + self.umap_bbknn = coords.get("bbknn") + self.umap_scanorama = coords.get("scanorama") + self.umap_seurat = coords.get("seurat") + + if subset_name is not None: + self.subset_name = subset_name + else: + self.subset_name = "Subset" + + sns.set_style("ticks") + + def df_get(self, subset, clustering, coords, category = None): + df = pd.DataFrame({ + "Subset" : np.repeat(subset, len(clustering)), + "UMAP 1" : coords[:, 0], + "UMAP 2" : coords[:, 1] + }) + subsets = natsorted(np.unique(subset)) + df["Subset"] = pd.Categorical( + np.repeat(subset, len(clustering)), categories=subsets, ordered=True + ) + df["Clustering"] = pd.Categorical( + clustering, categories=category, ordered=True + ) + return df + + def umap_df(self): + subset_list = [ + "bbknn", + "harmony", + "scanorama", + "scvi", + "seurat" + ] + clustering_list = [ + self.clustering_bbknn, + self.clustering_harmony, + self.clustering_scanorama, + self.clustering_scvi, + self.clustering_seurat + ] + clustering_unique = natsorted(np.unique(np.concatenate(clustering_list))) + coords_list = [ + self.umap_bbknn, + self.umap_harmony, + self.umap_scanorama, + self.umap_scvi, + self.umap_seurat + ] + umap_dfs = [ + self.df_get(i, j, k, category = clustering_unique) for i, j, k in zip( + subset_list, + clustering_list, + coords_list + ) + ] + self.umap_concat = pd.concat(umap_dfs) + + def umap_plot(self, show_plot = False): + self.umap_df() + palette = cc.glasbey_bw[0:len(np.unique(self.umap_concat["Clustering"]))] + self.umap_plt = sns.FacetGrid( + self.umap_concat, + col = "Subset", + col_wrap = 3, + hue = "Clustering", + palette = palette + ) + self.umap_plt.map( + sns.scatterplot, + "UMAP 1", + "UMAP 2", + s = 5, + alpha = 0.5 + ) + self.umap_plt.add_legend(markerscale = 3, title = self.subset_name) + if show_plot is True: + return self.umap_plt + + def save_umap(self, save_dir, dpi = 300): + try: + self.umap_plt.savefig( + save_dir, + dpi = dpi + ) + except: + self.umap_plot() + self.umap_plt.savefig( + save_dir, + dpi = dpi + ) \ No newline at end of file