diff --git a/.RData b/.RData new file mode 100644 index 0000000..fff6b6d Binary files /dev/null and b/.RData differ diff --git a/.Rhistory b/.Rhistory new file mode 100644 index 0000000..764522e --- /dev/null +++ b/.Rhistory @@ -0,0 +1,5 @@ +setwd("/home/lklossok/Morin_Lab/GAMBLR.predict") +devtools::load_all() +devtools::document() +devtools::load_all() +q() diff --git a/DESCRIPTION b/DESCRIPTION index b2c3098..77e5dc7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,7 +15,7 @@ Description: Collection of functions and helpers to classify different B-cell ly License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.2 Depends: R (>= 3.5.0) LazyData: true diff --git a/NAMESPACE b/NAMESPACE index 93e5d72..9cdcf5b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,22 +1,55 @@ # Generated by roxygen2: do not edit by hand +export(DLBCLone_KNN) +export(DLBCLone_KNN_predict) +export(DLBCLone_load_optimized) export(DLBCLone_optimize_params) +export(DLBCLone_predict_mixture_model) +export(DLBCLone_save_optimized) +export(DLBCLone_summarize_model) +export(DLBCLone_train_mixture_model) export(DLBCLone_train_test_plot) +export(assemble_genetic_features) +export(basic_umap_scatterplot) export(classify_bl) export(classify_dlbcl) export(classify_fl) export(complete_missing_from_matrix) export(construct_reduced_winning_version) +export(make_alluvial) export(make_and_annotate_umap) +export(make_neighborhood_plot) +export(make_umap_scatterplot) export(massage_matrix_for_clustering) -export(predict_single_sample) +export(optimize_outgroup) +export(optimize_purity) +export(predict_single_sample_DLBCLone) +export(prepare_single_sample_DLBCLone) +export(process_votes) +export(report_accuracy) +export(stacked_bar_plot) +export(summarize_all_ssm_status) export(tabulate_ssm_status) export(weighted_knn_predict_with_conf) +import(ComplexHeatmap) +import(FNN) import(GAMBLR.data) import(GAMBLR.helpers) +import(caret) +import(circlize) import(dplyr) +import(ggExtra) +import(ggalluvial) +import(ggplot2) +import(ggrepel) +import(ggside) +import(grid) +import(mclust) +import(purrr) import(randomForest, except = c("combine")) import(readr) +import(rlang) +import(stringr) import(tibble) import(tidyr) import(tidyselect) diff --git a/R/dlbclass.R b/R/dlbclass.R index 38ee54d..1be5596 100644 --- a/R/dlbclass.R +++ b/R/dlbclass.R @@ -28,13 +28,16 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB stop("Please provide either a mutations file or a mutation_data data frame.") } }else{ - mutation_data = read.table(mutations,sep="\t",row.names = 1,header=1) + mutation_data = read.table(mutations_file,sep="\t",row.names = 1,header=1) } # Transpose data if 'MYD88' is in row names if ("MYD88" %in% rownames(mutation_data)) { - mutation_data <- t(mutation_data) %>% as.data.frame() %>% column_to_rownames("sample_id") + mutation_data <- t(mutation_data) %>% as.data.frame() + #rownames(mutation_data) <- NULL + #print(head(mutation_data)) + #mutation_data = mutation_data %>% column_to_rownames("sample_id") } mutation_data = mutation_data %>% dplyr::select(-any_of(c("PLOIDY","PURITY","COO")),-ends_with("CCF")) #print(colnames(mutation_data)) @@ -72,34 +75,72 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB #mutation_data <- mutation_data[, !colnames(mutation_data) %in% genes_to_drop] # Aggregate specific features into vectors - print(grep("BCL6",colnames(mutation_data),value=T)) + + genes = c("BCL6","BCL2") BCL6_ALT <- rowSums(select(mutation_data, any_of(c("BCL6_SV","SV.BCL6", "BCL6"))), na.rm = TRUE) - NOTCH2_vec <- rowSums(select(mutation_data, any_of(c("NOTCH2","NOTCH2HOTSPOT", "SPEN", "DTX1"))), na.rm = TRUE) - M88O_vec <- rowSums(select(mutation_data, any_of(c("MYD88","MYD88.OTHER", "TNFAIP3", "TNIP1", "BCL10", "NFKBIE"))), na.rm = TRUE) + NOTCH2_genes <- c("NOTCH2","NOTCH2HOTSPOT", "SPEN", "DTX1") + genes = c(genes,NOTCH2_genes) + NOTCH2_vec <- rowSums(select(mutation_data, any_of(NOTCH2_genes)), na.rm = TRUE) - CD70_vec <- rowSums(select(mutation_data, any_of(c("CD70", "FAS", "CD58", "B2M", "FADD", "HLA.B","HLA-B"))), na.rm = TRUE) + MYD88_genes = c("MYD88","MYD88.OTHER", "TNFAIP3", "TNIP1", "BCL10", "NFKBIE") + genes = c(genes,MYD88_genes) + M88O_vec <- rowSums(select(mutation_data, any_of(MYD88_genes)), na.rm = TRUE) + CD70_genes = c("CD70", "FAS", "CD58", "B2M", "FADD", "HLA.B","HLA-B") + CD70_vec <- rowSums(select(mutation_data, any_of(CD70_genes)), na.rm = TRUE) + genes = c(genes,CD70_genes) # Additional vectors for other clusters #BCL2_combined <- rowSums(mutation_data[, c("BCL2", "SV.BCL2")], na.rm = TRUE) BCL2_combined <- rowSums(select(mutation_data, (starts_with("BCL2"))), na.rm = TRUE) + HIST_genes = c("HIST1H2AC","HIST1H2ACHOTSPOT", + "HIST1H1E", "HIST1H1EHOTSPOT","HIST1H1B", + "HIST1H2AM","HIST1H2AMHOTSPOT", "HIST1H1C", + "HIST1H1CHOTSPOT","HIST1H1D", "HIST1H1DHOTSPOT", + "HIST1H2BC","HIST1H2BCHOTSPOT") + genes = c(genes,HIST_genes) + SGK1_genes <- c("SGK1", "TET2", "SGK1HOTSPOT","NFKBIA","NFKBIAHOTSPOT", + "STAT3","STAT3HOTSPOT", "PTPN6", "BRAF","BRAFHOTSPOT", "KRAS", "KRASHOTSPOT", + "CD83","CD83HOTSPOT", "SF3B1","SF3B1HOTSPOT", "CD274", "MEF2C","MEF2CHOTSPOT", "KLHL6","KLHL6HOTSPOT", "CXCR4","CXCR4HOTSPOT", "PTEN", + "RAC2", "SESN3", "SOCS1","SOCS1HOTSPOT", "METAP1D") + DUSP2_genes <- c("DUSP2", "DUSP2HOTSPOT", "ZFP36L1","ZFP36L1HOTSPOT", "CRIP1", + "ACTB", "LTB", "YY1", "PABPC1","PABPC1HOTSPOT") - #CREBBP_vec <- rowSums(mutation_data[, c("CREBBP", "EZH2", "KMT2D", "EP300")], na.rm = TRUE) + TBL1XR1_genes <- c("TBL1XR1", "PIM1", "PIM1HOTSPOT", + "PRDM1","PRDM1HOTSPOT", + "ETV6","ETV6HOTSPOT", "ZC3H12A", + "BTG1", "BTG1HOTSPOT", "BTG2", "BTG2HOTSPOT", + "IGLL5", "IGLL5HOTSPOT", + "TMSB4X","TMSB4XHOTSPOT", + "GRHPR","GRHPRHOTSPOT","HLA.C","HLA-C","HLA-CHOTSPOT", + "MYD88", "TOX", "LYN", + "POU2F2","POU2F2HOTSPOT", + "IKZF3","IKZF3HOTSPOT", + "HLA.A", "ZFP36L1","HLA-A","HLA-AHOTSPOT", + "CARD11", "CARD11HOTSPOT","SF3B1", + "HLA.B","HLA-B","HLA-BHOTSPOT", "IRF2BP2", "OSBPL10", "ATP2A2", "PIM2", "IRF4", "BCL11A", + "METAP1D", "ETS1", "CCDC27") + CREBBP_genes <- c("CREBBP", "EZH2", "KMT2D", "EP300") + genes = c(genes,SGK1_genes,DUSP2_genes,TBL1XR1_genes,CREBBP_genes) CREBBP_vec <- rowSums(select(mutation_data, starts_with("CREBBP"), starts_with("EZH2"), starts_with("KMT2D"), starts_with("EP300")), na.rm = TRUE) - + PTEN_genes = c("PTEN") + C1_genes = c("UBE2A", "TMEM30A", "ZEB2", "GNAI2", "X5P.AMP", "POU2F2", "IKZF3", + "EBF1", "LYN", "BCL7A", "CXCR4", "CCDC27", "TUBGCP5", "SMG7", "RHOA", "BTG2") + TP53_genes = c("TP53") + GNA13_genes = c("GNA13", "TNFRSF14", "MAP2K1", "MEF2B", "IRF8", "HVCN1", + "GNAI2", "MEF2C", "SOCS1", "EEF1A1", "RAC2", + "POU2AF1") + genes = c(genes,PTEN_genes,C1_genes,TP53_genes,GNA13_genes) if(include_cn){ - C1_vec4 <- rowSums(select(mutation_data, any_of(c("UBE2A", "TMEM30A", "ZEB2", "GNAI2", "X5P.AMP", "POU2F2", "IKZF3", "X3Q28.DEL", - "EBF1", "LYN", "BCL7A", "CXCR4", "CCDC27", "TUBGCP5", "SMG7", "RHOA", "BTG2"))), na.rm = TRUE) - TP53_biallelic <- rowSums(mutation_data[, c("TP53", "X17P.DEL")], na.rm = TRUE) + C1_vec4 <- rowSums(select(mutation_data, any_of(C1_genes)), na.rm = TRUE) + TP53_biallelic <- rowSums(mutation_data[, TP53_genes,drop=FALSE], na.rm = TRUE) X21Q_AMP <- mutation_data[, "X21Q.AMP"] - GNA13_vec <- rowSums(mutation_data[, c("GNA13", "TNFRSF14", "MAP2K1", "MEF2B", "IRF8", "HVCN1", - "GNAI2", "MEF2C", "SOCS1", "EEF1A1", "RAC2", "X12Q.AMP", - "POU2AF1", "X6Q14.1.DEL")], na.rm = TRUE) + GNA13_vec <- rowSums(mutation_data[, GNA13_genes], na.rm = TRUE) Sum_C2_ARM <- rowSums(mutation_data[, c("X17P.DEL", "X21Q.AMP", "X11Q.AMP", "X6P.AMP", "X11P.AMP", "X6Q.DEL", "X7P.AMP", "X13Q.AMP", "X7Q.AMP", "X3Q.AMP", "X5P.AMP", "X18P.AMP", "X3P.AMP", "X19Q.AMP", "X9Q.AMP", "X12P.AMP", "X12Q.AMP")], na.rm = TRUE) @@ -109,7 +150,7 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB "X18Q23.DEL", "X19P13.3.DEL", "X13Q34.DEL", "X7Q22.1.AMP", "X10Q23.31.DEL", "X9P24.1.AMP", "X3Q28.AMP", "X11Q23.3.AMP", "X17Q24.3.AMP", "X3Q28.DEL", "X13Q14.2.DEL", "X18Q21.32.AMP", "X19Q13.32.DEL", "X6P21.1.AMP", "X18Q22.2.AMP", "EP300", "ZNF423", "CD274")], na.rm = TRUE) - PTEN <- rowSums(mutation_data[, c("PTEN", "X10Q23.31.DEL", "X13Q14.2.DEL")], na.rm = TRUE) + PTEN <- rowSums(mutation_data[, PTEN_genes, drop=FALSE], na.rm = TRUE) Sum_C5_CNA <- rowSums(mutation_data[, c("X18Q.AMP", "X3Q.AMP", "X3P.AMP", "X19Q13.42.AMP", "X6Q21.DEL", "X18P.AMP", "X19Q.AMP", "X8Q12.1.DEL", "X6Q14.1.DEL", "X19P13.2.DEL", "X9P21.3.DEL", "X18Q21.32.AMP", "X18Q22.2.AMP", "X1Q42.12.DEL", "X1Q32.1.AMP", "X6P21.33.DEL")], na.rm = TRUE) @@ -131,7 +172,7 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB starts_with("BTG2")), na.rm = TRUE) TP53_biallelic <- rowSums(mutation_data[, c("TP53"),drop=FALSE], na.rm = TRUE) - PTEN <- rowSums(mutation_data[, c("PTEN"),drop=FALSE], na.rm = TRUE) + PTEN <- rowSums(select(mutation_data,any_of(PTEN_genes))) #GNA13_vec <- rowSums(mutation_data[, c("GNA13", "TNFRSF14", "MAP2K1", "MEF2B", "IRF8", "HVCN1", # "GNAI2", "MEF2C", "SOCS1", "EEF1A1", "RAC2", # "POU2AF1")], na.rm = TRUE) @@ -149,40 +190,12 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB SV_MYC <- select(mutation_data, any_of(c("SV.MYC", "MYC","MYC_SV"))) %>% rowSums(na.rm = TRUE) - Hist_comp <- rowSums(select(mutation_data, any_of(c("HIST1H2AC","HIST1H2ACHOTSPOT", - "HIST1H1E", "HIST1H1EHOTSPOT","HIST1H1B", - "HIST1H2AM","HIST1H2AMHOTSPOT", "HIST1H1C", - "HIST1H1CHOTSPOT","HIST1H1D", "HIST1H1DHOTSPOT", - "HIST1H2BC","HIST1H2BCHOTSPOT"))), na.rm = TRUE) - SGK1_vec <- rowSums(select(mutation_data, any_of(c("SGK1", "TET2", "SGK1HOTSPOT","NFKBIA","NFKBIAHOTSPOT", - "STAT3","STAT3HOTSPOT", "PTPN6", "BRAF","BRAFHOTSPOT", "KRAS", "KRASHOTSPOT", - "CD83","CD83HOTSPOT", "SF3B1","SF3B1HOTSPOT", "CD274", "MEF2C","MEF2CHOTSPOT", "KLHL6","KLHL6HOTSPOT", "CXCR4","CXCR4HOTSPOT", "PTEN", - "RAC2", "SESN3", "SOCS1","SOCS1HOTSPOT", "METAP1D"))), na.rm = TRUE) - #SGK1_vec <- rowSums(mutation_data[, c("SGK1", "TET2", "NFKBIA", "STAT3", "PTPN6", "BRAF", "KRAS", - # "CD83", "SF3B1", "CD274", "MEF2C", "KLHL6", "CXCR4", "PTEN", - # "RAC2", "SESN3", "SOCS1", "METAP1D")], na.rm = TRUE) - #DUSP2_vec <- rowSums(mutation_data[, c("DUSP2", "ZFP36L1", "CRIP1", "ACTB", "LTB", "YY1", "PABPC1")], na.rm = TRUE) - DUSP2_vec <- rowSums(select(mutation_data, any_of(c("DUSP2", "DUSP2HOTSPOT", "ZFP36L1","ZFP36L1HOTSPOT", "CRIP1", "ACTB", "LTB", "YY1", "PABPC1","PABPC1HOTSPOT"))), na.rm = TRUE) + Hist_comp <- rowSums(select(mutation_data, any_of(HIST_genes)), na.rm = TRUE) + SGK1_vec <- rowSums(select(mutation_data, any_of(SGK1_genes)), na.rm = TRUE) + DUSP2_vec <- rowSums(select(mutation_data, any_of(DUSP2_genes)), na.rm = TRUE) + - #TBL1XR1_vec <- rowSums(mutation_data[, c("TBL1XR1", "PIM1", "PRDM1", "ETV6", "ZC3H12A", "BTG1", "BTG2", - # "IGLL5", "TMSB4X", "GRHPR", "HLA.C", "MYD88", "TOX", "LYN", - # "POU2F2", "IKZF3", "HLA.A", "ZFP36L1", "CARD11", "SF3B1", - # "HLA.B", "IRF2BP2", "OSBPL10", "ATP2A2", "PIM2", "IRF4", "BCL11A", - # "METAP1D", "ETS1", "CCDC27")], na.rm = TRUE) - TBL1XR1_vec <- rowSums(select(mutation_data, any_of(c("TBL1XR1", "PIM1", "PIM1HOTSPOT", - "PRDM1","PRDM1HOTSPOT", - "ETV6","ETV6HOTSPOT", "ZC3H12A", - "BTG1", "BTG1HOTSPOT", "BTG2", "BTG2HOTSPOT", - "IGLL5", "IGLL5HOTSPOT", - "TMSB4X","TMSB4XHOTSPOT", - "GRHPR","GRHPRHOTSPOT","HLA.C","HLA-C","HLA-CHOTSPOT", - "MYD88", "TOX", "LYN", - "POU2F2","POU2F2HOTSPOT", - "IKZF3","IKZF3HOTSPOT", - "HLA.A", "ZFP36L1","HLA-A","HLA-AHOTSPOT", - "CARD11", "CARD11HOTSPOT","SF3B1", - "HLA.B","HLA-B","HLA-BHOTSPOT", "IRF2BP2", "OSBPL10", "ATP2A2", "PIM2", "IRF4", "BCL11A", - "METAP1D", "ETS1", "CCDC27"))), na.rm = TRUE) + TBL1XR1_vec <- rowSums(select(mutation_data, any_of(TBL1XR1_genes)), na.rm = TRUE) MYD88_L265P_CD79B <- rowSums(select(mutation_data, any_of(c("MYD88.L265P","MYD88HOTSPOT", "CD79B","CD79BHOTSPOT"))), na.rm = TRUE) @@ -220,9 +233,6 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB C1_vec4 = C1_vec4, CD70_vec = CD70_vec, TP53_biallelic = TP53_biallelic, - #X21Q_AMP = X21Q_AMP, - #Sum_C2_ARM = Sum_C2_ARM, - #Sum_C2_FOCAL = Sum_C2_FOCAL, BCL2_combined = BCL2_combined, CREBBP_vec = CREBBP_vec, GNA13_vec = GNA13_vec, @@ -231,10 +241,8 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB Hist_comp = Hist_comp, SGK1_vec = SGK1_vec, DUSP2_vec = DUSP2_vec, - #CN_2P16_1_AMP = CN_2P16_1_AMP, TBL1XR1_vec = TBL1XR1_vec, MYD88_L265P_CD79B = MYD88_L265P_CD79B - #Sum_C5_CNA = Sum_C5_CNA ) } @@ -252,5 +260,10 @@ construct_reduced_winning_version <- function(mutations_file = "inst/extdata/DLB return(list(reduced=reduced_data, full=full_data, no_cn=data_no_cn, - feature_names=list(ssm_features=ssm,hotspot_features=hotspot,cnv_features=cnv,sv_features=sv))) + ssm_genes = genes, + feature_names=list(ssm_features=ssm, + hotspot_features=hotspot, + cnv_features=cnv, + sv_features=sv + ))) } diff --git a/R/plotting.R b/R/plotting.R new file mode 100644 index 0000000..4fd8980 --- /dev/null +++ b/R/plotting.R @@ -0,0 +1,1304 @@ +#' Heatmap visualization of mutations in nearest neighbors for a sample +#' +#' Generates a heatmap of feature values for the nearest neighbors of a specified sample, +#' based on a DLBCLone model object. This visualization helps to inspect the feature profiles +#' of samples most similar to the query sample. +#' +#' @param this_sample_id Character. The sample ID for which to plot the nearest neighbor heatmap. +#' @param DLBCLone_model A DLBCLone model object, which can be the output of \code{DLBCLone_optimize_params}, \code{DLBCLone_KNN}, or \code{predict_single_sample_DLBCLone}. +#' @param truth_column The column name in the predictions data frame that contains the ground truth labels (default: "lymphgen"). +#' @param metadata_cols Optional character vector of additional metadata columns to include in the heatmap annotations. +#' @param clustering_distance The distance metric to use for clustering rows. Default is "binary". +#' @param font_size Font size for heatmap text (default: 14). +#' +#' @return A ComplexHeatmap object showing the feature matrix for the nearest neighbors of the sample. +#' +#' @details +#' - For models of type \code{DLBCLone_optimize_params}, uses the \code{neighbors} and \code{lyseq_status} fields. +#' - For models of type \code{DLBCLone_KNN}, uses the \code{unlabeled_neighbors} field. +#' - The function extracts the feature matrix rows corresponding to the nearest neighbors of the specified sample, +#' and plots a heatmap of features with nonzero values. +#' +#' @import dplyr tidyr ComplexHeatmap grid circlize tibble +#' +#' @examples +#' # Assuming 'model' is a DLBCLone model and 'sample_id' is a valid sample: +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' nearest_neighbor_heatmap( +#' this_sample_id = "CCS_0680_lst", +#' DLBCLone_model = predict_single, +#' font_size = 10 +#' ) +#' } +#' +nearest_neighbor_heatmap <- function( + this_sample_id, + DLBCLone_model, + truth_column = "lymphgen", + metadata_cols = NULL, + clustering_distance = "binary", + font_size = 14 +){ + + if(!missing(DLBCLone_model) && "type" %in% names(DLBCLone_model)){ + if(DLBCLone_model$type == "DLBCLone_optimize_params"){ + neighbor_df = DLBCLone_model$neighbors + lyseq_status = DLBCLone_model$lyseq_status + }else if(DLBCLone_model$type == "DLBCLone_KNN"){ + + if(!"unlabeled_neighbors" %in% names(DLBCLone_model)){ + print(names(DLBCLone_model)) + stop("DLBCLone_model must be the output of DLBCLone_KNN with predict_unlabeled = TRUE") + }else if(is.null(DLBCLone_model$unlabeled_neighbors)){ + #no neighbors found for any of the incoming samples + message("No neighbors found for any sample. Returning NULL.") + return(NULL) + } + neighbor_df = DLBCLone_model$unlabeled_neighbors + neighbor_transpose = filter(neighbor_df,sample_id==this_sample_id) + if(nrow(neighbor_transpose) == 0){ + message("No neighbors found for sample ", this_sample_id, ". Returning NULL.") + return(NULL) + } + neighbor_transpose = neighbor_transpose %>% t() + #deal with fewer than K neighbours + neighbor_transpose = neighbor_transpose[!is.na(neighbor_transpose)] + pred_column = "DLBCLone_ko" + + }else if(DLBCLone_model$type == "predict_single_sample_DLBCLone"){ + DLBCLone_model$predictions <- DLBCLone_model$predictions %>% + filter(sample_id %in% this_sample_id) + + neighbor_ids = DLBCLone_model$predictions %>% + pull(neighbor_id) %>% + strsplit(",") %>% + unlist() + + neighbor_df <- DLBCLone_model$features_df %>% + rownames_to_column("sample_id") %>% + filter(sample_id %in% neighbor_ids) %>% + column_to_rownames("sample_id") + + neighbor_transpose = DLBCLone_model$features_df %>% + filter(rownames(DLBCLone_model$features_df) == this_sample_id) + + if(nrow(neighbor_transpose) == 0){ + message("No features found for sample ", this_sample_id, ". Returning NULL.") + return(NULL) + } + + neighbor_transpose = neighbor_transpose %>% t() + #deal with fewer than K neighbours + neighbor_transpose = neighbor_transpose[!is.na(neighbor_transpose)] + + pred_column = "predicted_label" + } + }else{ + stop("DLBCLone_model must be the output of DLBCLone_optimize_params, DLBCLone_KNN or predict_single_sample_DLBCLone") + } + + # Gather the base columns + cols_to_select <- c("sample_id", truth_column) + + # Add any non-null metadata columns + extra_cols <- c(pred_column,metadata_cols) + extra_cols <- extra_cols[!is.null(extra_cols)] + cols_to_select <- c(cols_to_select, extra_cols) + + xx = DLBCLone_model$features_df[neighbor_transpose,] + + if(any(is.na(rownames(xx)))){ + print(neighbor_transpose) + stop("something went wrong. Some samples are missing from features_df") + } + + top = max(xx) + mid = top/2 + col_fun = circlize::colorRamp2(c(0, mid, top), c("white", "#FFB3B3", "red")) + + if(DLBCLone_model$type == "predict_single_sample_DLBCLone"){ + + row_df = DLBCLone_model$optimized_predictions %>% + select(all_of(cols_to_select)) %>% + filter(sample_id %in% rownames(neighbor_df)) + + if(!this_sample_id %in% row_df$sample_id){ + new_row <- DLBCLone_model$predictions %>% + filter(sample_id %in% this_sample_id) %>% + mutate( + !!sym(truth_column) := NA, + !!sym(pred_column) := DLBCLone_model$predictions$predicted_label + ) + + # ensure extra metadata columns exist in the new row + for(col in extra_cols){ + if(!col %in% names(new_row)){ + new_row[[col]] <- NA + } + } + + new_row <- new_row %>% + select(all_of(cols_to_select)) + row_df <- bind_rows(row_df, new_row) + + }else{ + # create an NA row with all needed columns + missing_row <- as.data.frame(setNames(rep(NA, length(cols_to_select)), cols_to_select)) + missing_row$sample_id <- this_sample_id + row_df <- bind_rows(row_df, missing_row) + } + + }else{ + + #row_df = DLBCLone_model$predictions %>% + row_df = left_join(DLBCLone_model$predictions, DLBCLone_model$metadata) %>% + select(all_of(cols_to_select)) %>% + filter(sample_id %in% rownames(xx)) + + if(!this_sample_id %in% row_df$sample_id){ + if("unlabeled_predictions" %in% names(DLBCLone_model)){ + print("=====") + row_df = bind_rows( + row_df, + select( + DLBCLone_model$unlabeled_predictions, + sample_id, + !!sym(truth_column), + !!sym(pred_column) + ) %>% + filter(sample_id %in% rownames(xx)) + ) + print(row_df) + sample_class = filter(DLBCLone_model$unlabeled_predictions, sample_id == this_sample_id) %>% + pull(!!sym(pred_column)) + print(sample_class) + + }else{ + sample_class = NULL + row_df = bind_rows( + row_df, + tibble( + sample_id = this_sample_id, + !!rlang::sym(truth_column) := NA_character_, + !!rlang::sym(pred_column) := NA_character_) + ) + } + } + } + + feats_df <- DLBCLone_model$features_df %>% + rownames_to_column("sample_id") %>% + filter(sample_id %in% row_df$sample_id) %>% + column_to_rownames("sample_id") + + row_df = row_df %>% + column_to_rownames("sample_id") + + anno_colours = get_gambl_colours() + + anno_list = list() + + anno_list[[truth_column]] <- anno_colours + + for(col in extra_cols){ + anno_list[[col]] <- anno_colours + } + + row_anno = rowAnnotation( + df = row_df[rownames(feats_df),,drop=FALSE], + col = anno_list, + annotation_name_gp = gpar(fontsize = font_size), + show_legend = TRUE + ) + + title_text = paste( + "Sample", + this_sample_id, + "classified as", + pred_column + ) + + ht <- Heatmap( + feats_df[,colSums(feats_df)>0], + col = col_fun, + column_names_gp = gpar(fontsize = font_size), + right_annotation = row_anno, + clustering_distance_rows = clustering_distance, + show_heatmap_legend = FALSE, + column_title = title_text + ) + + draw( + ht, + heatmap_legend_side = "bottom", + annotation_legend_side = "bottom" + ) +} + +#' Basic UMAP Scatterplot +#' +#' Generates a simple UMAP scatterplot for visualizing sample clustering or separation. +#' +#' @param optimized Data frame containing at least V1, V2, sample_id, and grouping columns. +#' @param plot_samples Optional character vector of sample_ids to label in the plot +#' @param colour_by Column name to color points by. Defaults to `truth_column`. +#' @param truth_column Name of the truth/ground-truth column (default: "lymphgen"). +#' @param pred_column Name of the predicted-class column (default: "DLBCLone_ko"). +#' @param other_label Label used for the outgroup/unclassified class (default: "Other"). +#' @param title Plot title. +#' @param use_plotly Logical; if FALSE and `plot_samples` provided, draw static labels. +#' @param custom_colours Optional named vector of colors for groups; falls back to `get_gambl_colours()`. +#' +#' @return A ggplot object. +#' +#' @import dplyr ggplot2 ggExtra +#' @export +#' +#' @examples +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) +#' +#' basic_umap_scatterplot( +#' make_umap$df, #the data frame containing V1 and V2 from UMAP +#' plot_samples = "some_sample_ID", +#' colour_by = "DLBCLone_ko") +#' } +#' +basic_umap_scatterplot <- function(optimized, + plot_samples = NULL, + colour_by = NULL, + truth_column = "lymphgen", + pred_column = "DLBCLone_ko", + other_label = "Other", + title = "UMAP based on selected features", + use_plotly = TRUE, + custom_colours = NULL) { + stopifnot(all(c("V1","V2","sample_id") %in% colnames(optimized))) + stopifnot(is.data.frame(optimized)) + stopifnot(all(c("V1", "V2") %in% colnames(optimized))) + message("colour_by: ", colour_by) + colour_by <- colour_by %||% truth_column + + # dynamic label text (e.g., "lymphgen" and "DLBCLone_ko") + optimized_label <- optimized %>% + mutate( + label = paste0( + sample_id, + "\n", truth_column, ":\n", .data[[truth_column]], + " ", pred_column, ":\n", .data[[pred_column]] + ) + ) + + # points to label + label_points <- dplyr::filter(optimized_label, sample_id %in% (plot_samples %||% character())) %>% + mutate(label_x = V1 + 0.75, label_y = V2 - 0.75) + + # base mapping: color by chosen column + aes_cols <- aes( + x = V1, y = V2, + sample_id = sample_id, + color = .data[[colour_by]] + ) + + # draw outgroup first (by truth column) to keep it visually underneath + p <- ggplot() + p <- p + + geom_point( + data = dplyr::filter(optimized, .data[[truth_column]] == other_label), + mapping = aes_cols + ) + + geom_point( + data = dplyr::filter(optimized, .data[[truth_column]] != other_label), + mapping = aes_cols + ) + + # palette + pal <- custom_colours %||% get_gambl_colours() + p <- p + scale_color_manual(values = pal) + + guides(color = guide_legend(title = if (identical(colour_by, truth_column)) "Original Class" else "Predicted Class")) + + # optional static labels when not using plotly + if (!is.null(plot_samples) && length(plot_samples) > 0 && !isTRUE(use_plotly)) { + p <- p + + geom_segment( + data = label_points, + aes(x = V1, y = V2, xend = label_x, yend = label_y), + arrow = arrow(length = unit(0.02, "npc")), + color = "black" + ) + + geom_label( + data = label_points, + aes(x = label_x, y = label_y, label = label), + size = 3, + fill = "white", + label.size = 0.25 + ) + } + + p <- p + + labs(title = title) + + theme_minimal() + + return(p) +} + + +#' Summarize and Export DLBCLone Model Results +#' +#' Generates and saves a set of summary plots and tables for a DLBCLone model, including UMAP scatterplots, alluvial plots, and oncoplots. +#' Results are saved as PDF files in a directory named after the provided base name. +#' +#' @param base_name Character. The base name (and directory) for saving output files. +#' @param optimized_model List. The output from DLBCLone optimization, containing predictions, features, and metadata. +#' +#' @details +#' - Creates a directory for results if it does not exist. +#' - Saves UMAP scatterplots for all samples and for non-"Other" samples. +#' - Generates alluvial plots for different DLBCLone predictions +#' - Exports an oncoplot summarizing mutation and classification results. +#' - Uses `make_umap_scatterplot`, `make_alluvial`, and `prettyOncoplot` for visualization. +#' +#' @return No return value. Side effect: writes multiple PDF files to disk. +#' +#' @import dplyr ggplot2 ComplexHeatmap rlang +#' @export +#' +#' @examples +#' \dontrun{ +#' DLBCLone_summarize_model("Full_geneset_unweighted", optimized_model) +#'} +#' +DLBCLone_summarize_model = function( + base_name, + optimized_model +){ + base_dir = here::here() + full_dir = paste0(base_dir,"/",base_name) + if(!dir.exists(full_dir)){ + dir.create(full_dir) + } + #save the predictions + pred_out = paste0(full_dir,"/DLBCLone_predictions_bulk.tsv") + print(paste("Saving predictions to",pred_out)) + to_save = optimized_model$predictions %>% + select(sample_id,!!sym(optimized_model$truth_column),DLBCLone_i:DLBCLone_wo,confidence:other_score,V1,V2) + + write_tsv(to_save, + file = pred_out) + umap1 = paste0(full_dir,"/UMAP_all.pdf") + cairo_pdf(umap1,width=8,height=8) + p = make_umap_scatterplot( + optimized_model$df, + colour_by = optimized_model$truth_column, + title="all samples, projected" + ) + print(p) + dev.off() + + umap2 = paste0(full_dir,"/UMAP_no_Other.pdf") + cairo_pdf(umap2,width=8,height=8) + p = make_umap_scatterplot( + optimized_model$df, + colour_by = optimized_model$truth_column, + drop_other = T, + title="non-Other samples, projected" + ) + print(p) + dev.off() + cairo_pdf(paste0(full_dir,"/Alluvial_DLBCLone_io.pdf"),width=8,height=12) + p = make_alluvial( + optimized_model, + pred_column = "DLBCLone_io", + truth_column = optimized_model$truth_column + ) + print(p) + dev.off() + cairo_pdf(paste0(full_dir,"/Alluvial_DLBCLone_wo.pdf"),width=8,height=12) + p = make_alluvial( + optimized_model, + pred_column = "DLBCLone_wo", + truth_column = optimized_model$truth_column + ) + print(p) + dev.off() + cairo_pdf(paste0(full_dir,"/Oncoplot.pdf"),width=16,height=14) + mc = c("lymphgen","DLBCLone_wo","DLBCLone_io","DLBCLone_w","DLBCLone_i") + sc = c("DLBCLone_wo","DLBCLone_io","lymphgen","DLBCLone_w","DLBCLone_i","confidence") + if("DLBClass" %in% colnames(optimized_model$predictions)){ + mc = c(mc,"DLBClass") + sc = c(sc,"DLBClass") + } + mc = mc[mc %in% colnames(optimized_model$predictions)] + sc = sc[sc %in% colnames(optimized_model$predictions)] + prettyOncoplot( + all_maf_with_s, + these_samples_metadata =optimized_model$predictions, + genes = colnames(optimized_model$features), + minMutationPercent = 1, + metadataColumns = mc, + sortByColumns = sc, + numericMetadataColumns = "confidence", + simplify_annotation = T, + cluster_rows=T + ) + dev.off() +} + +#' @title Make Neighborhood Plot +#' @description +#' Generates a UMAP plot highlighting the neighborhood of a given sample, showing its nearest neighbors and their group assignments. +#' +#' @param single_sample_prediction_output A list containing prediction results and annotation data frames. +#' Must include elements \code{prediction} (data frame with prediction results) and \code{anno_df} (data frame with UMAP coordinates and annotations). +#' @param training_predictions The equivalent data frame of prediction results for all training samples (e.g. optimized_model$df) +#' @param this_sample_id Character. The sample ID for which the neighborhood plot will be generated. +#' @param prediction_in_title Logical. If \code{TRUE}, includes the predicted label in the plot title. +#' @param add_circle Plot will include a circle surrounding the set of neighbors. Set to FALSE to disable. +#' @param label_column Does nothing, i.e. this is not currently working. +#' +#' @return A \code{ggplot2} object representing the UMAP plot with the selected sample and its neighbors highlighted. +#' +#' @details +#' The function extracts the nearest neighbors of the specified sample, draws segments connecting the sample to its neighbors, and colors points by group (e.g., lymphgen subtype). The plot title can optionally include the predicted label. +#' +#' @import dplyr ggplot2 rlang ggExtra +#' @export +#' +#' @examples +#' +#' # Assuming 'optimization_result' is the output of DLBCLone_optimize_params +#' # and 'output' is the result of DLBCLone_predict_single_sample +#' # on sample_id "SAMPLE123": +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' predict_single <- predict_single_sample_DLBCLone( +#' test_df = optimize_params$features[1,], +#' train_metadata = dlbcl_meta, +#' optimized_model = optimize_params +#' ) +#' +#' make_neighborhood_plot( +#' single_sample_prediction_output = predict_single, +#' training_predictions = optimize_params$df, +#' this_sample_id = "SAMPLE123", +#' prediction_in_title = TRUE, +#' add_circle = TRUE +#' ) +#' } +#' +make_neighborhood_plot <- function( + single_sample_prediction_output, + training_predictions, + this_sample_id, + prediction_in_title = TRUE, + add_circle = TRUE, + label_column = "DLBCLone_io", + point_size = 0.5, + point_alpha = 0.9, + line_alpha = 0.9 +){ + + circleFun <- function(center = c(0,0),diameter = 1, npoints = 100){ + r = diameter / 2 + tt <- seq(0,2*pi,length.out = npoints) + xx <- center[1] + r * cos(tt) + yy <- center[2] + r * sin(tt) + return(data.frame(x = xx, y = yy)) + } + if(missing(training_predictions)){ + #training_predictions = single_sample_prediction_output$anno_df + training_predictions = + left_join( + select(single_sample_prediction_output$anno_df,sample_id,lymphgen), + select(single_sample_prediction_output$df,sample_id,V1,V2) + ) + }else if(missing(single_sample_prediction_output)){ + #Just plot the single sample in the context of the rest based on the optimization + single_sample_prediction_output = list() + single_sample_prediction_output[["predictions"]] = filter(training_predictions, sample_id==this_sample_id) + single_sample_prediction_output[["anno_df"]] = training_predictions + }else{ + single_sample_prediction_output$predictions = filter(single_sample_prediction_output$predictions, sample_id==this_sample_id) + } + xmin = min(training_predictions$V1, na.rm = TRUE) + xmax = max(training_predictions$V1, na.rm = TRUE) + ymin = min(training_predictions$V2, na.rm = TRUE) + ymax = max(training_predictions$V2, na.rm = TRUE) + #extract the sample_id for all the nearest neighbors with non-Other labels + my_neighbours = filter(single_sample_prediction_output$predictions,sample_id == this_sample_id) %>% + pull(neighbor_id) %>% strsplit(.,",") %>% unlist() + + #set up links connecting each neighbor to the sample's point + links_df = filter(training_predictions,sample_id %in% my_neighbours) %>% mutate(group=lymphgen) + my_x = filter(single_sample_prediction_output$anno_df,sample_id==this_sample_id) %>% pull(V1) + my_y = filter(single_sample_prediction_output$anno_df,sample_id==this_sample_id) %>% pull(V2) + if(prediction_in_title){ + title = paste(this_sample_id,pull(single_sample_prediction_output$predictions,!!sym(label_column))) + if(single_sample_prediction_output$predictions[[label_column]] == "Other" && single_sample_prediction_output$predictions$predicted_label !="Other"){ + title = paste(title,"(",single_sample_prediction_output$predictions$predicted_label,")") + } + + }else{ + title = this_sample_id + } + links_df = mutate(links_df,my_x=my_x,my_y=my_y) + links_df = links_df %>% select(V1,V2,my_x,my_y,group) %>% mutate(length = sqrt((V1 - my_x)^2 + (V2 - my_y)^2)) # Euclidean distance + + pp=ggplot(mutate(training_predictions,group=lymphgen),aes(x=V1,y=V2,colour=group)) + + geom_point(alpha=point_alpha,size=point_size) + + geom_segment(data=links_df,aes(x=V1,y=V2,xend=my_x,yend=my_y),alpha=line_alpha) + + scale_colour_manual(values=get_gambl_colours()) + + ggtitle(title) + + xlim(c(xmin,xmax)) + + ylim(c(ymin,ymax)) + + theme_minimal() + if(add_circle){ + #add a circle around the sample + d = max(links_df$length)*2.1 # adding a 10% spacer + circle = circleFun(c(my_x,my_y),diameter=d,npoints=100) + pp = pp + geom_path(data=circle,aes(x=x,y=y),colour="black",alpha=1,inherit.aes=FALSE) + } + return(pp) +} + +#' Make UMAP scatterplot +#' +#' @param df Data frame containing the UMAP coordinates and annotations. +#' @param truth_column Name of the column containing the labels to be plotted (default: "lymphgen"). +#' @param drop_composite If TRUE: removes composite labels from the lymphgen column. +#' @param drop_other If TRUE: removes "Other" and "NOS" labels from the lymphgen column. +#' @param high_confidence If TRUE: filters the data to include only samples with confidence > 0.7. +#' @param custom_colours Custom color palette for the plot. If not provided, uses default GAMBL colors. +#' @param add_labels If TRUE: adds labels to the points based on the median coordinates of each group. +#' @param title Title for the plot. Default: NULL. +#' +#' @returns A ggplot object representing the UMAP scatterplot with marginal histograms. +#' +#' @import dplyr ggplot2 ggExtra ggside rlang +#' @export +#' +#' @examples +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) +#' +#' make_umap_scatterplot( +#' make_umap$df, +#' drop_other = F +#' ) +#' } +#' +make_umap_scatterplot = function( + df, + truth_column = "lymphgen", + drop_composite = TRUE, + drop_other = FALSE, + high_confidence = FALSE, + custom_colours, + add_labels = FALSE, + title = NULL +){ + xmin = min(df$V1,na.rm=TRUE) + xmax = max(df$V1,na.rm=TRUE) + ymin = min(df$V2,na.rm=TRUE) + ymax = max(df$V2,na.rm=TRUE) + if(!missing(custom_colours)){ + cols = custom_colours + }else{ + cols = get_gambl_colours() + } + if(drop_composite){ + df = filter(df,!is.na(!!sym(truth_column)),!grepl("COMP",!!sym(truth_column))) + } + if(drop_other){ + df = filter(df,!is.na(!!sym(truth_column)),!!sym(truth_column)!="Other",!!sym(truth_column)!="NOS") + } + if(high_confidence){ + df = filter(df,Confidence > 0.7) + } + if(add_labels){ + labels = group_by(df,!!sym(truth_column)) %>% + summarise(median_x = median(V1),median_y = median(V2)) + } + unique_lg = unique(df[[truth_column]]) + if(any(!unique_lg %in% names(cols))){ + missing = unique_lg[!unique_lg %in% names(cols)] + print(paste("missing colour for:",paste(missing,collapse=","))) + } + p = ggplot(df,aes(x=V1,y=V2,colour=!!sym(truth_column),label=cohort)) + + geom_point(alpha=0.8) + + scale_colour_manual(values=cols) + theme_Morons() + + guides(colour = guide_legend(nrow = 1)) + + xlim(xmin,xmax) + + ylim(ymin,ymax) + + if(!is.null(title)){ + p = p + ggtitle(title) + } + if(add_labels){ + p = p + geom_label_repel(data=labels,aes(x=median_x,y=median_y,label=!!sym(truth_column))) + } + ggMarginal(p,groupColour = TRUE,groupFill=TRUE) +} + +#' Report Classification Accuracy and Per-Class Metrics +#' +#' Computes overall accuracy, balanced accuracy, and sensitivity for predicted vs. true class labels. +#' Optionally excludes samples assigned to the "Other" class from accuracy calculations. +#' +#' @param predictions Data frame containing predicted and true class labels. +#' @param truth Name of the column with true class labels (default: "lymphgen"). +#' @param pred Name of the column with predicted class labels (default: "predicted_label"). +#' @param per_group Logical; if TRUE, computes per-group accuracy metrics. +#' @param metric Character; type of accuracy to report ("accuracy" supported). +#' @param verbose Whether to print verbose outputs to console +#' +#' @return A list with: +#' \item{no_other}{Accuracy excluding samples assigned to "Other"} +#' \item{per_class}{Balanced accuracy per class} +#' \item{per_class_sensitivity}{Sensitivity per class} +#' \item{overall}{Overall accuracy including all samples} +#' +#' @details +#' - Uses confusion matrices to compute accuracy metrics. +#' - Excludes "Other" class for no_other accuracy. +#' - Returns per-class metrics for further analysis. +#' +#' @import dplyr caret rlang +#' @export +#' +#' @examples +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' result <- report_accuracy(predictions_df) +#' result$overall +#' result$per_class +#' } +#' +report_accuracy = function( + predictions, + truth="lymphgen", + pred="DLBCLone_io", + per_group = FALSE, + metric = "accuracy", + verbose = FALSE +){ + if("BN2" %in% predictions[[pred]]){ + truth = "lymphgen" + } + all_classes = unique(predictions[[truth]]) + no_other_pred = filter(predictions,lymphgen != "Other") + if(verbose){ + print(paste(nrow(no_other_pred),"non-Other samples were assigned a label")) + print(paste(filter(no_other_pred,!!sym(pred) !="Other") %>% nrow(),"non-Other samples were assigned to a non-Other class")) + } + + conf_matrix_no <- confusionMatrix( + factor(no_other_pred[[pred]],levels=all_classes), + factor(no_other_pred[[truth]],levels=all_classes) + ) + #print(conf_matrix) + overall = conf_matrix_no$overall[["Accuracy"]] + + if(metric == "accuracy"){ + #bal_acc_no <- conf_matrix$byClass[, "Balanced Accuracy"] + acc_no = overall + }else{ + stop("unsupported metric") + } + + conf_matrix <- confusionMatrix( + factor(predictions[[pred]],levels=unique(predictions[[truth]])), + factor(predictions[[truth]],levels=unique(predictions[[truth]])) + ) + + if(metric == "accuracy"){ + bal_acc <- conf_matrix$byClass[, "Balanced Accuracy"] + sensitivity = conf_matrix$byClass[, "Sensitivity"] + }else{ + stop("unsupported metric") + } + overall = conf_matrix$overall[["Accuracy"]] + + return(list( + no_other=acc_no, + per_class=bal_acc, + mean_balanced_accuracy = mean(bal_acc,na.rm=TRUE), + per_class_sensitivity = sensitivity, + overall=overall, + confusion_matrix_no_other=conf_matrix_no, + confusion_matrix=conf_matrix + )) +} + +#' Create an Alluvial Plot Comparing Original and Predicted Classifications +#' +#' This function generates a detailed alluvial plot to visualize the concordance and discordance between original (e.g., Lymphgen) and predicted (e.g., DLBCLone) class assignments for samples. +#' It supports annotation of concordance rates, per-group accuracy, unclassified rates, and flexible labeling and coloring options. +#' +#' @param optimized List containing prediction results and metadata, typically output from a DLBCLone optimization function. +#' @param count_excluded_as_other Logical; if TRUE, samples excluded due to missing features are counted as "Other" in the plot. +#' @param title Plot title (default: empty string). +#' @param group_order Character vector specifying the order of groups/classes for axes and coloring. +#' @param add_accuracy_to_title Logical; if TRUE, adds accuracy/concordance rate to the plot title. +#' @param accuracy_per_group Logical; if TRUE, computes and displays per-group accuracy. +#' @param accuracy_type Type of accuracy to report (default: "sensitivity"). +#' @param truth_name Name for the original class column (default: "Lymphgen"). +#' @param truth_column Name for the original class column in the predictions data frame (default: "lymphgen"). +#' @param pred_name Name for the predicted class column (default: "DLBCLone"). +#' @param pred_column Name of the column in predictions to use for the predicted class (default: "predicted_label_optimized"). +#' @param nudge Amount to nudge labels horizontally (default: 0.03). +#' @param box_nudge Amount to nudge label boxes (default: 0.15). +#' @param min_flip_n Minimum number of samples for a flow to be labeled (default: 10). +#' @param add_percent Logical; if TRUE, adds percent concordance to labels. +#' @param add_unclass_rate Logical; if TRUE, adds unclassified rate annotation to the plot. +#' @param denom Denominator for stratum/flow width scaling (default: 20). +#' @param label_size Size of label text (default: 3). +#' @param label_lock_y Logical; if TRUE, locks label movement to the x-axis only. +#' @param label_line_flip_colour Logical; controls color assignment for label lines. +#' @param label_box_flip_colour Logical; controls color assignment for label boxes. +#' @param concordant_label_relative_pos Position for concordant labels: 0 (left), 0.5 (middle), or 1 (right). +#' @param verbose Whether to print verbose outputs to console +#' +#' @return A ggplot2 object representing the alluvial plot. +#' +#' @details +#' - Visualizes flows between original and predicted classes, highlighting concordant and discordant assignments. +#' - Annotates concordance rate, per-group accuracy, and unclassified rate as specified. +#' - Supports flexible labeling, coloring, and axis ordering for publication-quality plots. +#' +#' @import dplyr ggplot2 ggalluvial rlang tidyr +#' @export +#' +#' @examples +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' make_alluvial(optimize_params) +#' } +#' +make_alluvial <- function( + optimized, + count_excluded_as_other = FALSE, + title = "", + group_order = NULL, + add_accuracy_to_title = TRUE, + accuracy_per_group = TRUE, + accuracy_type = "sensitivity", + truth_name = "Lymphgen", + truth_column = "lymphgen", + pred_name = "DLBCLone", + pred_column = "predicted_label_optimized", + nudge = 0.03, + box_nudge = 0.15, + min_flip_n = 10, + add_percent = TRUE, + add_unclass_rate = TRUE, + denom = 20, + label_size=3, + label_lock_y=TRUE, + label_line_flip_colour=TRUE, + label_box_flip_colour=TRUE, + concordant_label_relative_pos = 0.5, #we only handle 0, 0.5 and 1 + verbose = FALSE +){ + predictions = optimized$predictions + if(is.null(group_order)){ + group_order = optimized$truth_classes + print("setting group order:") + } + + #print(group_order) + + if(is.null(truth_column) && !is.null(optimized$truth_column)){ + truth_column = optimized$truth_column + } + if (accuracy_per_group) { + accuracies <- report_accuracy( + predictions, + truth = truth_column, + pred = pred_column, + per_group = accuracy_per_group + ) + } + #print("Done accuracyies") + if(count_excluded_as_other){ + excluded_meta = optimized$sample_metadata_no_features %>% + mutate(!!pred_column := "Other") + predictions = bind_rows(excluded_meta,predictions) + } + if("total_samples_available" %in% names(optimized)){ + full_denominator = optimized$total_samples_available + }else{ + full_denominator = nrow(predictions) + } + + xx <- predictions %>% + #group_by(lymphgen, predicted_label_optimized) %>% + #summarize(num = n(), .groups = "drop") %>% + rename( + !!truth_name := !!sym(truth_column), + #!!new_name := predicted_label_optimized + !!pred_name := !!sym(pred_column) + ) + #print("Done renaming") + xx <- xx %>% + group_by(!!sym(truth_name), !!sym(pred_name)) %>% + summarize(num = n(), .groups = "drop") + + + grid <- expand_grid( + !!sym(truth_name) := group_order, + !!sym(pred_name) := group_order + ) + #print(colnames(grid)) + #print(colnames(xx)) + xx <- grid %>% + left_join(xx, by = c(truth_name, pred_name)) %>% + mutate(num = replace_na(num, 0)) + #print("Done joining") + xx[[truth_name]] <- factor(xx[[truth_name]], levels = group_order) + xx[[pred_name]] <- factor(xx[[pred_name]], levels = group_order) + + xx = xx %>% mutate(pair = paste(!!sym(truth_name),!!sym(pred_name),sep="-")) + xy = filter(xx,!!sym(truth_name)!="Other") %>% + mutate(match=ifelse(!!sym(truth_name)==!!sym(pred_name),TRUE,FALSE)) %>% + group_by(match) %>% + summarise(concordant=sum(num)) %>% + ungroup() %>% + mutate(total=sum(concordant)) %>% + mutate(percent=100*concordant/total) + + pc_conc = filter(xy,match==TRUE) %>% pull(percent) %>% round(.,1) + title = paste0(title," Concordance (non-Other): ",pc_conc,"%") + bacc = round(accuracies$mean_balanced_accuracy,4) + + # One side sort + + xx <- xx %>% + group_by(!!sym(truth_name)) %>% + arrange(!!sym(pred_name), .by_group = TRUE) %>% + ungroup() %>% + mutate(flow_id = row_number(), flow_label = ifelse(num > min_flip_n, as.character(num), "")) + + xx = xx %>% + group_by(!!sym(pred_name)) %>% + arrange(!!sym(truth_name),.by_group=TRUE) %>% + ungroup() %>% + mutate(flow_id2 = row_number()) + + lodes <- ggalluvial::to_lodes_form( + xx, + key = "axis", + axes = c(truth_name, pred_name), + id = "flow_id" + ) %>% + mutate( + axis = as.numeric(axis), + nudge_dir = 0.1 + ) %>% + left_join(xx %>% select(flow_id, !!sym(truth_name), !!sym(pred_name)), by = "flow_id") + + stratum_bases <- lodes %>% + mutate(stratum = factor(stratum, levels = group_order)) %>% + group_by(axis, stratum) %>% + summarize(total = sum(num), .groups = "drop") %>% + arrange(axis, desc(stratum)) %>% + group_by(axis) %>% + mutate(stratum_base = cumsum(total) - total) %>% + ungroup() + + left_lodes <- lodes %>% + left_join(stratum_bases, by = c("axis", "stratum")) %>% + separate(pair,into=c("left_group","right_group"),sep = "-") %>% + mutate(right_group = factor(right_group,levels=levels(stratum))) %>% + mutate(left_group = factor(left_group,levels=levels(stratum))) %>% + group_by(axis, stratum) %>% + arrange( + desc(stratum), + desc(right_group), + num, + .by_group = TRUE + ) %>% + mutate( + flow_bottom = cumsum(num) - num, + flow_y = stratum_base + flow_bottom + num / 2 + ) %>% + ungroup() + + right_lodes <- lodes %>% + left_join(stratum_bases, by = c("axis", "stratum")) %>% + separate(pair,into=c("left_group","right_group"),sep = "-") %>% + mutate(right_group = factor(right_group,levels=levels(stratum))) %>% + mutate(left_group = factor(left_group,levels=levels(stratum))) %>% + group_by(axis, stratum) %>% + arrange( + desc(right_group), + desc(left_group), + num, + .by_group = TRUE) %>% + mutate( + flow_bottom = cumsum(num) - num, + flow_y = stratum_base + flow_bottom + num / 2 + ) %>% + ungroup() + + #print("HERE") + lodes_mean_left = group_by(left_lodes,right_group,left_group) %>% + mutate(y=mean(flow_y)) %>% ungroup() %>% + arrange(right_group,left_group) %>% + mutate(nudge_dir = -nudge) + #print("HERE") + lodes_mean_right = group_by(right_lodes,left_group,right_group) %>% + mutate(y=mean(flow_y)) %>% ungroup() %>% + arrange(left_group,right_group) %>% + mutate(nudge_dir = nudge) + + if(concordant_label_relative_pos == 0){ + lodes_mean = lodes_mean_left %>% mutate(y=flow_y) %>% + filter(axis==1) + #position on far left + }else if(concordant_label_relative_pos == 1){ + lodes_mean = lodes_mean_right %>% + mutate(y=flow_y) %>% + filter(axis==2) + #position on far right + }else if(concordant_label_relative_pos == 0.5){ + lodes_mean = bind_rows(filter(right_lodes,axis==2),filter(left_lodes,axis==1)) %>% + group_by(left_group,right_group) %>% + mutate(y=mean(flow_y)) %>% ungroup() %>% + arrange(left_group,right_group) %>% + filter(axis==2) %>% + mutate(nudge_dir = 0) + #middle + }else{ + stop("unhandled value for concordant_label_relative_pos. Provide 0 (left), 0.5 (middle) or 1 (right)") + } + #print("HERE") + lodes_right = filter(lodes_mean_right,axis==2,left_group!=right_group) %>% + arrange(flow_id2) %>% + mutate(nudge_dir = -nudge) %>% + mutate(left_char=as.character(left_group),right_char = as.character(right_group)) + + + lodes_left = filter(lodes_mean_left,axis==1,left_group!=right_group) %>% + arrange(flow_id) %>% + mutate(left_char=as.character(left_group),right_char = as.character(right_group)) %>% + ungroup() + + lodes_left = lodes_left %>% + mutate(segment_colour = if (!label_line_flip_colour) left_char else right_char) %>% + mutate(label_fill = if(label_box_flip_colour) right_group else left_group) %>% + mutate(nudge_dir = nudge) + + lodes_right = lodes_right %>% + mutate(segment_colour = if (label_line_flip_colour) left_char else right_char) %>% + mutate(label_fill = if(!label_box_flip_colour) right_group else left_group) %>% + mutate(nudge_dir = -nudge) + + lodes_match = filter(lodes_mean,left_group==right_group) %>% + mutate(nudge_dir=0) + #print(lodes_match) + #print("HERE!") + + xx_denom <- predictions %>% + group_by(!!sym(truth_column)) %>% + summarize(denom = n(), .groups = "drop") %>% + rename( + stratum = !!sym(truth_column) + ) + #print("HERE!!") + + if(add_percent){ + lodes_match = left_join(lodes_match,xx_denom,by="stratum") %>% + mutate(percent = round(100*num/denom,1)) %>% + mutate(flow_label = paste0(flow_label," (",percent,"%",")")) + } + if(add_unclass_rate){ + right_other = filter(lodes_mean_right,right_group=="Other",axis==2) %>% + summarize(total=sum(num)) %>% + mutate(unclass= round(100 * total / full_denominator,1)) %>% + mutate(label=paste0(total,"\n (",unclass,"%)")) %>% + mutate(!!sym(truth_name) := "Other",!!sym(pred_name) := "Other") + #print(right_other) + unclass = pull(right_other,unclass) + title = paste0(title,", Classification rate: ",100-unclass,"%") + title = paste0(title," Bal Acc: ",bacc) + left_other = filter(lodes_mean_left,left_group=="Other",axis==2) %>% + summarize(total=sum(num)) %>% + mutate(unclass= round(100 * total / full_denominator,1)) %>% + mutate(label=paste0(total,"\n (",unclass,"%)")) %>% + mutate(!!sym(truth_name) := "Other",!!sym(pred_name) := "Other") + } + + pp = ggplot(xx, aes( + axis1 = !!sym(truth_name), + axis2 = !!sym(pred_name), + y = num + )) + + geom_alluvium(aes(fill = !!sym(truth_name)), width = 1 / denom) + + geom_stratum(width = 1 / denom, aes(fill = after_stat(stratum)), color = "black") + + geom_label( + data = filter(lodes_match,left_group==right_group,flow_label!=""), + stat = "identity", + inherit.aes = FALSE, + colour="white", + size=label_size, + aes( + x = concordant_label_relative_pos + 1, + y, + label = flow_label, + fill = right_group + ) + ) + + geom_label_repel( + data = filter(lodes_right,flow_label!=""), + min.segment.length = 0, + stat = "identity", + inherit.aes = FALSE, + direction = ifelse(label_lock_y,"x","both"), + colour="white", + size=label_size, + nudge_x = box_nudge, + aes( + x = axis + nudge_dir, + flow_y, + label = flow_label, + fill = label_fill, + segment.colour = segment_colour + ) + ) + + geom_label_repel( + data = filter(lodes_left,flow_label!=""), + min.segment.length = 0, + stat = "identity", + inherit.aes = FALSE, + direction = ifelse(label_lock_y,"x","both"), + colour="white", + size=label_size, + nudge_x = -box_nudge, + aes( + x = axis + nudge_dir, + flow_y, + label = flow_label, + fill = label_fill, + segment.colour = segment_colour + ) + ) + + scale_x_discrete(limits = c(truth_name, pred_name), expand = c(.1, .05)) + + scale_fill_manual(values = get_gambl_colours()) + + scale_colour_manual(values = get_gambl_colours(),aesthetics = c("color", "segment.color")) + + labs( + y = "Number of Samples", + x = NULL + ) + + theme_minimal() + + theme(legend.position = "none") + + ggtitle(title) + + if(add_unclass_rate){ + pp = pp + + geom_label_repel( + data=left_other, + direction = "x", + x=1, + y=20, + aes(fill=!!sym(truth_name),label=label) + ) + + geom_label_repel( + data=right_other, + direction = "x", + x=2, + y=20, + aes(fill=!!sym(pred_name),label=label) + ) + } + pp +} + +#' Create a stacked bar plot of top features per subtype +#' +#' This function generates a stacked bar plot showing the top features (genes) for each subtype based on their prevalence in the dataset. +#' +#' @param DLBCLone_model A DLBCLone model object, which can be the output of \code{DLBCLone_optimize_params}, or \code{DLBCLone_KNN} +#' @param truth_column Name of the column containing the true class labels (default: "lymphgen"). +#' @param truth_classes Vector of class labels to consider (default: c("BN2","EZB","MCD","ST2")). +#' @param method Method to determine top features: "common" for most common features, "chi_square" for subtype vs rest significance (default : "common"). +#' @param num_feats Number of top features to display per subtype (default: 10). +#' @param title Title for the plot (default: NULL). +#' +#' @return A ggplot2 object representing the stacked bar plot. +#' +#' @import dplyr ggplot2 tidyr rlang +#' @export +#' +#' @examples +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' stacked_bar_plot( +#' optimize_params, +#' method = "chi_square", +#' num_feats = 10, +#' title = "LymphGen" +#' ) +#' } +#' +stacked_bar_plot <- function( + DLBCLone_model, + truth_column = "lymphgen", + truth_classes = c("BN2","EZB","MCD","ST2"), + method = "common", + num_feats = 10, + title = NULL +){ + if(!missing(DLBCLone_model) && "type" %in% names(DLBCLone_model) && DLBCLone_model$type == "predict_single_sample_DLBCLone"){ + stop("DLBCLone_model must be the output of DLBCLone_optimize_params, or DLBCLone_KNN") + } + + if ("type" %in% names(DLBCLone_model) && DLBCLone_model$type == "DLBCLone_KNN") { + DLBCLone_model$features <- DLBCLone_model$features_df + } + + bad_cols <- colSums(DLBCLone_model$features) <= 0.02 * nrow(DLBCLone_model$features) + DLBCLone_model$features <- DLBCLone_model$features[, !bad_cols] + + annotated_feats <- DLBCLone_model$features %>% + rownames_to_column("sample_id") %>% + left_join( + DLBCLone_model$df %>% select(sample_id, !!sym(truth_column)), + by = "sample_id" + ) %>% + column_to_rownames("sample_id") + + annotated_feats[[truth_column]] <- as.factor(annotated_feats[[truth_column]]) + + gene_cols <- setdiff(colnames(annotated_feats), c("sample_id", truth_column)) + subtypes <- truth_classes + + top_genes_per_subtype <- list() + + if(method == "common"){ + + for(subtype in subtypes){ + subtype_samples <- annotated_feats[annotated_feats[[truth_column]] == subtype, ] + + gene_counts <- colSums(subtype_samples[, gene_cols, drop = FALSE]) + gene_ranking <- order(gene_counts, decreasing = TRUE) + + # Top genes for this subtype + top_genes_per_subtype[[subtype]] <- gene_cols[gene_ranking[1:num_feats]] + } + + }else if(method == "chi_square"){ + + for(subtype in subtypes){ + # Create binary class: subtype vs rest + y_binary <- factor(ifelse(annotated_feats[[truth_column]] == subtype, subtype, paste0("not_", subtype))) + + chi_stats <- numeric(length(gene_cols)) + + for(i in seq_along(gene_cols)){ + gene <- gene_cols[i] + tbl <- table(annotated_feats[[gene]], y_binary) + test <- suppressWarnings(chisq.test(tbl)) + chi_stats[i] <- test$statistic + } + + gene_ranking <- order(chi_stats, decreasing = TRUE) + + # Top genes for this subtype + top_genes_per_subtype[[subtype]] <- gene_cols[gene_ranking[1:num_feats]] + } + + }else{ + stop( + paste( + "Must select a valid method:", + "'common' - for most common features", + "'chi_square' - for subtype vs rest significance", + sep = "\n" + ) + ) + } + + plot_df <- bind_rows(lapply(names(top_genes_per_subtype), function(subtype){ + genes <- top_genes_per_subtype[[subtype]] + subtype_samples <- annotated_feats[annotated_feats[[truth_column]] == subtype, ] + n_subtype <- nrow(subtype_samples) # total samples in this subtype + + counts <- colSums(subtype_samples[, genes, drop = FALSE] > 0) # ensure 0/1 + plot_df <- tibble( + subtype = subtype, + gene = names(counts), + count = as.numeric(counts), + prop_gene = as.numeric(counts) / n_subtype + ) + })) + + plot_df <- plot_df %>% + group_by(subtype) %>% + mutate( + prop = count / sum(count), + rank = rank(-count, ties.method = "first") # rank genes within subtype + ) %>% + ungroup() + + # cumulative position for placing labels + plot_df <- plot_df %>% + group_by(subtype) %>% + arrange(rank) %>% + mutate( + cum_prop = cumsum(prop), + pos = cum_prop - prop / 2 + ) %>% + ungroup() + + colours <- get_gambl_colours() + n_colours <- max(plot_df$rank) + + fill_map <- plot_df %>% + group_by(subtype) %>% + mutate( + base_col = ifelse(!is.na(colours[subtype]), colours[subtype], "grey"), + ramp = list(colorRampPalette(c(first(base_col), "white"))(max(rank))), + fill_color = ramp[[1]][rank] + ) %>% + ungroup() + + # merge colors back into plot_df + plot_df$fill_color <- fill_map$fill_color + + ggplot(plot_df, aes(x = subtype, y = prop, fill = fill_color)) + + geom_bar(stat = "identity", color = "black", position = position_stack(reverse = TRUE)) + + geom_text( + aes(label = paste0(gene, " ", scales::percent(prop_gene, accuracy = 1))), + position = position_stack(vjust = 0.5, reverse = TRUE), + size = 3 + ) + + scale_fill_identity() + + labs( + title = paste(title, " Top", num_feats, "genes per subtype (method:", method, ")"), + x = "Subtype", + y = "Proportion of Samples with Mutation" + ) + + theme_minimal() + + theme( + axis.text.x = element_text(angle = 0, hjust = 1), + legend.position = "none" + ) +} diff --git a/R/umap.R b/R/umap.R index e35114b..03c6575 100644 --- a/R/umap.R +++ b/R/umap.R @@ -1,89 +1,461 @@ -#' Predict class for a single sample without using umap_transform +#' Summarize SSM (Somatic Single Nucleotide Mutation) Status Across Samples #' -#' @param test_df Data frame containing the mutation status of the test sample -#' @param train_df Data frame containing the mutation status of the training samples -#' @param best_params Data frame from DLBCLone_optimize_params with the best parameters -#' @param train_metadata Metadata for training samples with truth labels in lymphgen column -#' @param truth_classes Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2") -#' @param drop_unlabeled_from_training Set to TRUE to drop unlabeled samples from the training data -#' @param make_plot Set to TRUE to plot the UMAP projection and predictions +#' This function summarizes the mutation status for a set of genes +#' across multiple samples, separating mutations by class for genes specified +#' by \code{separate_by_class_genes} and counting the number of hits +#' in each sample per mutation category. #' -#' @returns a list of data frames with the predictions and the UMAP input +#' @param maf_df A data frame containing mutation annotation format (MAF) data, +#' with at least the following columns: +#' \code{Hugo_Symbol}, \code{Variant_Classification}, and \code{Tumor_Sample_Barcode}. +#' @param these_samples_metadata A data frame containing metadata for the samples, +#' with at least a \code{sample_id} column. +#' Any sample that does not have a matching sample_id in these_samples_metadata will be dropped. +#' @param genes_of_interest A character vector of gene symbols to include +#' in the summary. If missing, defaults to all Tier 1 B-cell lymphoma genes. +#' @param synon_genes (Optional) A character vector of gene symbols for which +#' synonymous mutations should be included. +#' @param silent_maf_df (Optional) A separate data frame containing silent mutation data if +#' the user doesn't want to pull silent mutation status from \code{maf_df}. +#' This argument is useful when you want to combine mutations from +#' the output of get_coding_ssm and get_ssm_by_region or get_ssm_by_gene +#' @param separate_by_class_genes (Optional) A character vector of +#' gene symbols for which mutations should be separated by class +#' (e.g., "Nonsense_Mutation", "Missense_Mutation"). +#' @param count_hits Logical; if \code{TRUE}, counts the number of mutations +#' per gene per sample. If \code{FALSE} (default), only presence/absence is recorded. +#' +#' @return A wide-format data frame (matrix) with samples as rows and +#' mutation types as columns. Each cell contains either the count of +#' mutations (if \code{count_hits = TRUE}) or a binary indicator (0/1) +#' for mutation presence. +#' +#' @details +#' - Mutations are grouped and optionally separated by mutation +#' class for genes specified in \code{separate_by_class_genes} +#' - Synonymous mutations can be counted as another separate feature for genes specified by \code{synon_genes} genes. +#' - The function simplifies mutation annotations and pivots the data to a wide format suitable for downstream analysis. +#' +#' @import dplyr tidyr tibble +#' @export +#' +#' @examples +#' # A basic example, using only the output of get_all_coding_ssm +#' # Since the only non-coding class this function handles is Silent, +#' # we will be missing most non-coding types such as Intron, UTR, Flank +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' sample_metadata = get_sample_metadata() %>% filter(seq_type!= "mrna") +#' +#' maf_data = get_all_coding_ssm(sample_metadata) +#' +#' mutation_matrix <- summarize_all_ssm_status( +#' maf_df = maf_data, +#' these_samples_metadata = sample_metadata, +#' genes_of_interest = c("TP53", "SGK1", "BCL2"), +#' synon_genes = c("BCL2"), +#' separate_by_class_genes = c("TP53","SGK1"), +#' count_hits = FALSE +#' ) +#' } +#' +summarize_all_ssm_status <- function(maf_df, + these_samples_metadata, + genes_of_interest, + synon_genes, + silent_maf_df, + separate_by_class_genes = NULL, + count_hits = FALSE){ + if(missing(genes_of_interest)){ + message("defaulting to all Tier 1 B-cell lymphoma genes") + genes_of_interest = filter(GAMBLR.data::lymphoma_genes, + DLBCL_Tier==1 | FL_Tier == 1 | BL_Tier == 1 ) %>% + pull(Gene) %>% unique() + } + + maf_df = filter(maf_df, + Hugo_Symbol %in% genes_of_interest) + if(missing(silent_maf_df)){ + silent_maf_df = maf_df + }else{ + silent_maf_df = filter(silent_maf_df,Hugo_Symbol %in% genes_of_interest) + } + if(!missing(these_samples_metadata)){ + maf_df = filter(maf_df, + Tumor_Sample_Barcode %in% these_samples_metadata$sample_id) + silent_maf_df = filter(silent_maf_df, + Tumor_Sample_Barcode %in% these_samples_metadata$sample_id) + } + if(!missing(synon_genes)){ + if(any(!synon_genes %in% silent_maf_df$Hugo_Symbol)){ + missing = synon_genes[!synon_genes %in% silent_maf_df$Hugo_Symbol] + not_missing = synon_genes[synon_genes %in% silent_maf_df$Hugo_Symbol] + + message("Warning: Some synonymous genes have no mutations in silent_maf_df") + message(paste(missing,collapse=", ")) + } + maf_nonsilent = filter(maf_df, Variant_Classification %in% vc_nonSynonymous) + maf_silent = filter(silent_maf_df, ! Variant_Classification %in% vc_nonSynonymous, Hugo_Symbol %in% synon_genes) + maf_df = bind_rows(maf_silent,maf_nonsilent) + } + #Simplify annotations + silent_types = c("Silent","Intron","5'UTR","3'UTR","5'Flank","3'Flank") + nonsense_types = c("Nonsense_Mutation","Frame_Shift_Del","Frame_Shift_Ins","Splice_Site") + missense_types = c("In_Frame_Del", + "In_Frame_Ins", + "Missense_Mutation", + "Splice_Region") + gene_mutations = mutate(maf_df, + mutation_type=case_when( + Variant_Classification %in% nonsense_types ~ "Nonsense_Mutation", + Variant_Classification %in% missense_types ~ "Missense_Mutation", + Variant_Classification %in% silent_types ~ "Silent", + TRUE ~ Variant_Classification) + ) %>% + mutate(mutation=paste(Hugo_Symbol,mutation_type,sep=":")) + separated_coding_maf = filter(gene_mutations, + Hugo_Symbol %in% genes_of_interest, + Hugo_Symbol %in% separate_by_class_genes, + mutation_type != "Silent") + + unseparated_coding_maf = filter(gene_mutations, + Hugo_Symbol %in% genes_of_interest, !Hugo_Symbol %in% separate_by_class_genes, mutation_type != "Silent") %>% + mutate(mutation = paste(Hugo_Symbol,"Coding",sep=":")) + separated_silent_maf = filter(gene_mutations,Hugo_Symbol %in% synon_genes, mutation_type == "Silent") %>% + mutate(mutation = paste(Hugo_Symbol,"Silent",sep=":")) + gene_mutations = bind_rows(separated_coding_maf, unseparated_coding_maf,separated_silent_maf) + + mutation_distinct = select(gene_mutations,mutation,Tumor_Sample_Barcode) %>% + mutate(mutated = 1) %>% + group_by(Tumor_Sample_Barcode,mutation,mutated) %>% + count() %>% ungroup() + if(count_hits){ + mutation_distinct = mutation_distinct %>% select(-mutated) %>% rename(mutated=n) %>% distinct() + } else{ + mutation_distinct = mutation_distinct %>% select(-n) %>% distinct() + } + mutation_wide = pivot_wider(mutation_distinct, names_from = "mutation",values_from = "mutated",values_fill = 0) %>% + column_to_rownames("Tumor_Sample_Barcode") + return(mutation_wide) +} + +#' Assemble genetic features for UMAP input +#' +#' This function assembles a matrix of genetic features for each sample, including mutation status, +#' aSHM counts, and structural variant status for BCL2, BCL6, and MYC. It supports both genome and capture sequencing types. +#' +#' @param these_samples_metadata Data frame with sample metadata, must include seq_type and sample_id. +#' @param metadata_columns Columns in metadata to use for SV status (default: c("bcl2_ba","bcl6_ba","myc_ba")). +#' @param genes Vector of gene symbols to include. +#' @param synon_genes Vector of gene symbols for synonymous mutations. +#' @param maf_with_synon MAF data frame including synonymous mutations. +#' @param hotspot_genes Vector of hotspot genes. +#' @param genome_build Genome build version (default: "grch37"). +#' @param sv_value Value to assign for SV presence (default: 3). +#' @param synon_value Value to assign for synonymous mutations (default: 1). +#' @param coding_value Value to assign for coding mutations (default: 2). +#' @param include_ashm Logical; if TRUE, includes aSHM counts in the output matrix (default: TRUE). +#' @param annotated_sv Data frame with annotated structural variants. +#' @param include_GAMBL_sv Logical; if TRUE, includes GAMBL structural variants in the output matrix (default: TRUE). +#' @param review_hotspots Logical; if TRUE, reviews hotspot genes (default: TRUE). +#' @param verbose Defaults to FALSE +#' +#' @return Matrix of assembled features for each sample. +#' +#' @import dplyr tibble tidyr #' @export +#' +#' @examples +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' ordered_genes = filter(lymphoma_genes,DLBCL_Tier==1 | FL_Tier==1|BL_Tier==1) %>% pull(Gene) +#' synon_genes = unique(grch37_ashm_regions$gene) +#' synon_genes = synon_genes[synon_genes %in% ordered_genes] +#' +#' all_sv_anno = get_combined_sv(these_samples_metadata = dlbcl_meta) #' -predict_single_sample = function(test_df, - train_df, - train_metadata, - best_params, - truth_classes = c("EZB","MCD","ST2","N1","BN2"), - drop_unlabeled_from_training=TRUE, - make_plot = TRUE){ - train_metadata_use = filter(train_metadata,lymphgen %in% truth_classes) - train_metadata_notuse = filter(train_metadata,!lymphgen %in% truth_classes) - - if(nrow(test_df)>1){ - message("Warning: you have supplied more than one sample to test with. Will proceed with all") - } - combined_mutation_status_df = bind_rows(test_df,train_df) - if(any(rownames(test_df) %in% train_metadata_use$sample_id)){ - stop("one or more samples overlap with your training data!") - } - placeholder_meta = data.frame(sample_id = rownames(test_df)) - train_metadata_use= bind_rows(placeholder_meta,train_metadata_use) - - outs = make_and_annotate_umap(df=combined_mutation_status_df, - min_dist = 0, - n_neighbors = 55, - init="random", - n_epochs = 1500, - seed=best_params$seed, - metadata=train_metadata_use, - ret_model=T, - metric="cosine", - join_column="sample_id", - na_vals = best_params$na_option) - - test_coords = dplyr::filter(outs$df, - #!sample_id %in% rownames(test_df), - is.na(lymphgen)) %>% select(sample_id,V1,V2) %>% +#' all_sv_anno = annotate_sv(all_sv_anno) +#' +#' all_full_status = assemble_genetic_features( +#' these_samples_metadata = dlbcl_meta, +#' metadata_columns = c("BCL2_SV","BCL6_SV"), +#' synon_genes = synon_genes, +#' synon_value = 1, +#' coding_value = 2, +#' genes = ordered_genes, +#' hotspot_genes = c("MYD88"), +#' maf_with_synon = all_maf_with_s, +#' include_ashm = F, +#' sv_value = 2, +#' verbose=F, +#' annotated_sv = all_sv_anno +#' ) +#' } +#' +assemble_genetic_features <- function(these_samples_metadata, + metadata_columns = c("bcl2_ba","bcl6_ba","myc_ba"), + genes, + synon_genes, + maf_with_synon, + hotspot_genes, + genome_build = "grch37", + sv_value = 3, + synon_value = 1, + coding_value = 2, + include_ashm = TRUE, + annotated_sv, + include_GAMBL_sv= TRUE, + review_hotspots = TRUE, + verbose = FALSE){ + + if(!"maf_data" %in% class(maf_with_synon)){ + warning(paste("maf_with_synon should be a maf_data object, but is not.", + "Proceeding anyway with genome_build = ", genome_build)) + }else{ + genome_build = attr(maf_with_synon,"genome_build") + } + if(include_ashm){ + stopifnot( genome_build == "grch37", "other genome builds are not yet supported for aSHM") + #TODO: consider shifting this function to GAMBLR.results + #TODO: ensure this supports both genome builds correctly. This is currently hard-coded + some_regions = GAMBLR.utils::create_bed_data( + GAMBLR.data::grch37_ashm_regions, + fix_names = "concat", + concat_cols = c("gene","region"), + sep="-") + #make the aSHM count matrix and combine if necessary + if ("genome" %in% these_samples_metadata$seq_type){ + ashm_matrix_genome <- get_ashm_count_matrix( + regions_bed = some_regions, + this_seq_type = "genome", + these_samples_metadata = these_samples_metadata + ) + + colnames(ashm_matrix_genome) = gsub("-.+","",colnames(ashm_matrix_genome)) + } + + if ("capture" %in% these_samples_metadata$seq_type){ + ashm_matrix_cap <- get_ashm_count_matrix( + regions_bed = some_regions, + this_seq_type = "capture", + these_samples_metadata = these_samples_metadata + ) + colnames(ashm_matrix_cap) = gsub("-.+","",colnames(ashm_matrix_cap)) + } + if ("genome" %in% these_samples_metadata$seq_type && "capture" %in% these_samples_metadata$seq_type){ + ashm_matrix = bind_rows(ashm_matrix_genome, ashm_matrix_cap) + }else if("genome" %in% these_samples_metadata$seq_type){ + ashm_matrix = ashm_matrix_genome + }else if("capture" %in% these_samples_metadata$seq_type){ + ashm_matrix = ashm_matrix_cap + + }else{ + stop("no eligible seq_type provided in these_samples_metadata") + } + } + + include_hotspots = ifelse(!missing(hotspot_genes), TRUE, FALSE) + + status_with_silent = get_coding_ssm_status( + these_samples_metadata = these_samples_metadata, + # drop all coding variants from this one + maf_data = maf_with_synon, + include_hotspots = include_hotspots, + genes_of_interest = hotspot_genes, + include_silent_genes = synon_genes[synon_genes %in% genes], + gene_symbols = genes + ) + status_with_silent = status_with_silent %>% column_to_rownames("sample_id") + + status_without_silent = get_coding_ssm_status( + these_samples_metadata = these_samples_metadata, + maf_data = maf_with_synon, + include_hotspots = include_hotspots, + genes_of_interest = hotspot_genes, + include_silent = FALSE, + gene_symbols = genes + ) + + status_without_silent = status_without_silent %>% column_to_rownames("sample_id") - pred = weighted_knn_predict_with_conf( - train_coords = dplyr::filter(outs$df, - #!sample_id %in% rownames(test_df), - !is.na(lymphgen)) %>% select(V1,V2), - train_labels = dplyr::filter(outs$df, - #!sample_id %in% rownames(test_df) - !is.na(lymphgen)) %>% pull(lymphgen), - test_coords = test_coords, - k=best_params$k, - conf_threshold = best_params$threshold, - na_label="Other", - use_weights = best_params$use_w, - ignore_top = FALSE) - pred$sample_id = rownames(test_coords) - - if(drop_unlabeled_from_training){ - pred = dplyr::filter(pred,sample_id %in% rownames(test_df)) - } - #print(head(outs$df)) - anno_out = left_join(pred,outs$df,by="sample_id") %>% - mutate(label = paste(sample_id,predicted_label,round(confidence,3))) - if(make_plot){ - pp = ggplot(outs$df,aes(x=V1,y=V2,colour=lymphgen,label=sample_id)) + - geom_point() + - geom_point(data = anno_out,aes(colour=predicted_label)) + - geom_label_repel(data = anno_out, - aes(x=V1,y=V2,label=label), - nudge_x=1,nudge_y=1,colour="black") + - scale_colour_manual(values=get_gambl_colours()) + - ggtitle(paste("k:",best_params$k,"seed:",best_params$seed)) - print(pp) - } - return(list(prediction = pred, umap_input = outs$features, model=outs$model)) + if (any(! colnames(status_with_silent) %in% colnames(status_without_silent))){ + print(colnames(status_with_silent)[!colnames(status_with_silent) %in% colnames(status_without_silent)]) + stop("some columns are missing from the status_without_silent matrix") + } + # Instead of just relying on the MAF(s) supplied by the user + if (include_ashm){ + ashm_matrix = select(ashm_matrix, any_of(colnames(status_with_silent))) %>% select(any_of(synon_genes)) + ashm_matrix[ashm_matrix>1]= 1 + + if(verbose){ + print(head(colSums(ashm_matrix))) + print(head(ashm_matrix[,c(1:10)])) + } + #fill in gaps from aSHM (other non-coding variants in the genes) + + missing = status_with_silent[rownames(ashm_matrix), + colnames(ashm_matrix)]==0 & + ashm_matrix[rownames(ashm_matrix), + colnames(ashm_matrix)] > 0 + + + fill = missing + fill[]=0 + fill[missing] = synon_value + status_with_silent[rownames(fill), + colnames(fill)] = fill + + } + + #ensure all columns in status_with_silent are present in status_without_silent + if (any(! colnames(status_without_silent) %in% colnames(status_with_silent))){ + print(colnames(status_without_silent)[!colnames(status_without_silent) %in% colnames(status_with_silent)]) + stop("some columns are missing from the status_with_silent matrix") + } + + # ensure column order is identical in the two matrices + + status_with_silent = status_with_silent[,colnames(status_without_silent)] + + status_combined = status_with_silent + status_without_silent + if(coding_value == 1){ + status_combined[status_combined > 1] = 1 + } + #TODO: generalize this to use the column names provided in metadata_columns + bcl2_id = these_samples_metadata[which(these_samples_metadata$bcl2_ba=="POS"),] %>% pull(sample_id) + bcl6_id = these_samples_metadata[which(these_samples_metadata$bcl6_ba=="POS"),] %>% pull(sample_id) + myc_id = these_samples_metadata[which(these_samples_metadata$myc_ba=="POS"),] %>% pull(sample_id) + + # include SV from provided annotated SV data frame + if(!missing(annotated_sv) & include_GAMBL_sv){ + annotated_sv = filter(annotated_sv,tumour_sample_id %in% these_samples_metadata$sample_id) + bcl2_sv_id = filter(annotated_sv,!is.na(partner),gene=="BCL2") %>% pull(tumour_sample_id) + bcl6_sv_id = filter(annotated_sv,!is.na(partner),gene=="BCL6") %>% pull(tumour_sample_id) + myc_sv_id = filter(annotated_sv,!is.na(partner),gene=="MYC") %>% pull(tumour_sample_id) + n_b = length(bcl2_id) + + bcl2_id = unique(c(bcl2_id,bcl2_sv_id)) + n_b = length(bcl2_id) + + bcl6_id = unique(c(bcl6_id,bcl6_sv_id)) + myc_id = unique(c(myc_id,myc_sv_id)) + } + + + status_combined[,"BCL2_SV"] = 0 + status_combined[bcl2_id,"BCL2_SV"] = sv_value + + status_combined[,"MYC_SV"] = 0 + status_combined[myc_id,"MYC_SV"] = sv_value + + status_combined[,"BCL6_SV"] = 0 + status_combined[bcl6_id,"BCL6_SV"] = sv_value + return(status_combined) + } -#' Plot the result of a lymphgen classification +#' Optimize the threshold for classifying samples as "Other" +#' +#' Performs a post-hoc evaluation of the classification of a sample as one of +#' the main classes vs the outgroup/unclassified label "Other" and returns the +#' optimal threshold for classifying a sample as "Other" based on the ground +#' truth provided in the true_labels vector. It evaluates the performance +#' of the classifier using a range of thresholds and returns the best threshold +#' based on the specified metric (balanced accuracy or accuracy). +#' This function is not generally meant to be called directly but rather is +#' a helper function used by DLBCLone_optimize_params. +#' +#' @param predicted_labels Vector of predicted labels for the samples +#' @param true_labels Vector of true labels for the samples +#' @param other_score Vector of scores for the "Other" class for each sample () +#' @param all_classes Vector of classes to use for training and testing. +#' Default: c("MCD","EZB","BN2","N1","ST2","Other") +#' @param maximize Metric to use for optimization. +#' Either "accuracy" (actual accuracy across all samples) or "balanced_accuracy" +#' (the mean of the balanced accuracy values across all classes). +#' Default: "balanced_accuracy" +#' @param exclude_other_for_accuracy Set to TRUE to exclude the +#' "Other" class from the 'lymphgen' column when calculating accuracy metrics +#' (passed to DLBCLone_optimize_params). Default: FALSE +#' @param verbose Set to TRUE to print additional information during the optimization process. +#' +#' @returns a list of data frames with the predictions and the UMAP input +#' +#' @import dplyr tidyr caret +#' @export +#' +optimize_outgroup <- function(predicted_labels, + true_labels, + other_score, + all_classes = c("MCD", + "EZB", + "BN2", + "N1", + "ST2", + "Other"), + maximize ="balanced_accuracy", + exclude_other_for_accuracy = FALSE, + cap_classification_rate = 1, + verbose = FALSE, + other_class = "Other"){ + rel_thresholds = seq(1,10,0.1) + sens_df = data.frame() + acc_df = data.frame() + predictions = data.frame(predicted_label=as.character(predicted_labels), + true_label=as.character(true_labels)) + + for(threshold in rel_thresholds){ + predictions_new = mutate(predictions, + predicted_label = ifelse(other_score < threshold, + predicted_label, + other_class)) + all_acc = report_accuracy(predictions_new, + truth="true_label", + pred="predicted_label", + drop_other=FALSE) + pred = factor(predictions_new[["predicted_label"]],levels=all_classes) + truth = factor(predictions_new[["true_label"]],levels=all_classes) + conf_matrix <- confusionMatrix(pred, truth) + + bal_acc <- conf_matrix$byClass[, "Balanced Accuracy"] + if(maximize == "balanced_accuracy"){ + bal_acc$average_accuracy = mean(bal_acc) + }else{ + bal_acc$average_accuracy = conf_matrix$overall[["Accuracy"]] + } + bal_acc$threshold = threshold + bal_acc$harmonic_mean = all_acc$harmonic_mean + bal_acc$classification_rate = all_acc$classification_rate + acc_df = bind_rows(acc_df,bal_acc) + sn <- conf_matrix$byClass[, "Sensitivity"] + sn$average_sensitivity = mean(sn) + sn$threshold = threshold + sn$harmonic_mean = all_acc$harmonic_mean + sn$classification_rate = all_acc$classification_rate + sens_df = bind_rows(sens_df,sn) + } + + + if(maximize %in% c("balanced_accuracy","accuracy")){ + acc_df = filter(acc_df, classification_rate <= cap_classification_rate) + best = slice_head(arrange(acc_df,desc(average_accuracy)),n=1) + }else if(maximize == "harmonic_mean"){ + acc_df = filter(sens_df, classification_rate <= cap_classification_rate) + best = slice_head(arrange(sens_df,desc(harmonic_mean)),n=1) + }else{ + best = slice_head(arrange(sens_df,desc(average_sensitivity)),n=1) + + } + + return(best) +} + +#' Plot the result of a DLBCLone classification #' #' @param test_df Data frame containing the test data with UMAP coordinates #' @param train_df Data frame containing the training data with UMAP coordinates @@ -98,10 +470,16 @@ predict_single_sample = function(test_df, #' @param title3 additional argument #' #' @returns a ggplot object +#' +#' @import ggplot2 dplyr ggrepel #' @export #' #' @examples #' #add the dataset name to the metadata if it's not already there (required for the plot work) +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' #' lymphgen_A53_DLBCLone$df$dataset = "GAMBL" #' #' DLBCLone_train_test_plot( @@ -112,6 +490,7 @@ predict_single_sample = function(test_df, #' details = lymphgen_A53_DLBCLone$best_params, #' classes = c("MCD","EZB","BN2","ST2","N1","A53","Other"), #' annotate_accuracy=TRUE,label_offset = 1) +#' } #' DLBCLone_train_test_plot = function(test_df, train_df, @@ -119,32 +498,38 @@ DLBCLone_train_test_plot = function(test_df, other_df, details, annotate_accuracy = FALSE, + classes = c("BN2","ST2","MCD","EZB","N1"), label_offset = 2, - title1="GAMBL", - title2="predicted_class_for_HighConf", - title3 ="predicted_class_for_Other"){ - - title = paste0("N_class:",details$num_classes," N_feats:",details$num_features," k=",details$k," threshold=",details$threshold," bacc=",round(details$accuracy,3)) - if("BN2" %in% classes){ - print(details) - acc_df = data.frame(lymphgen = c("N1","BN2","EZB","MCD","ST2","Other","A53"), - accuracy = c(details$N1_bacc, - details$BN2_bacc, - details$EZB_bacc, - details$MCD_bacc, - details$ST2_bacc, - details$Other_bacc, - details$A53_bacc)) - }else if("C1" %in% classes){ - acc_df = data.frame(lymphgen = c("C1","C2","C3","C4","C5"), - accuracy = c(details$C1_bacc, - details$C2_bacc, - details$C3_bacc, - details$C4_bacc, - details$C5_bacc)) - }else{ - stop("no labels to add?") + title1="Original Class", + title2="DLBCLone Predicted Class", + title3 ="DLBCLone Predicted Class (Other)",base_size = 1){ + title = "" + if(!missing(details)){ + title = paste0("N_class:",details$num_classes," N_feats:",details$num_features," k=",details$k," threshold=",details$threshold," bacc=",round(details$accuracy,3)) + + } + if(annotate_accuracy){ + if("BN2" %in% classes){ + acc_df = data.frame(lymphgen = classes, + accuracy = c( + details$BN2_bacc, + details$EZB_bacc, + details$MCD_bacc, + details$ST2_bacc, + details$Other_bacc, + details$A53_bacc)) + }else if("C1" %in% classes){ + acc_df = data.frame(lymphgen = c("C1","C2","C3","C4","C5"), + accuracy = c(details$C1_bacc, + details$C2_bacc, + details$C3_bacc, + details$C4_bacc, + details$C5_bacc)) + }else{ + stop("no labels to add?") + } + } # Add the predicted labels for Other (unclassified) cases, if provided if(!missing(other_df)){ @@ -153,18 +538,20 @@ DLBCLone_train_test_plot = function(test_df, mutate(predictions_df,dataset=title2,lymphgen=predicted_label), mutate(other_df,dataset=title3,lymphgen=predicted_label) ) + in_df = mutate(in_df,dataset = factor(dataset,levels=unique(c(unique(train_df$dataset),title1,title2,title3)))) }else{ in_df = bind_rows(train_df, test_df, mutate(predictions_df,dataset=title2,lymphgen=predicted_label) ) + in_df = mutate(in_df,dataset = factor(dataset,levels=unique(c(unique(train_df$dataset),title1,title2)))) } pp = ggplot(in_df) + geom_point(aes(x=V1,y=V2,colour=lymphgen),alpha=0.8) + scale_colour_manual(values=get_gambl_colours()) + facet_wrap(~dataset,ncol=1) + - theme_Morons() + ggtitle(title) + theme_Morons(base_size=base_size) + ggtitle(title) if(annotate_accuracy){ #add labels and set nudge direction based on what quadrant each group sits in centroids = filter(predictions_df,predicted_label %in% classes) %>% @@ -172,6 +559,7 @@ DLBCLone_train_test_plot = function(test_df, summarise(mean_V1=median(V1),mean_V2=median(V2)) %>% mutate(nudge_x=sign(mean_V1),nudge_y = sign(mean_V2)) %>% mutate(lymphgen=predicted_label) + #print(centroids) centroids = left_join(centroids,acc_df) %>% mutate(label=paste(lymphgen,":",round(accuracy,3))) @@ -185,10 +573,9 @@ DLBCLone_train_test_plot = function(test_df, geom_label_repel(data=filter(centroids,nudge_y > 0, nudge_x > 0), aes(x=mean_V1,y=mean_V2,label=label),fill="white",size=5,nudge_y = 1 * label_offset , nudge_x = 1 * label_offset) } - pp + pp + guides(colour = guide_legend(nrow = 1)) } - #' Run UMAP and attach result to metadata #' #' @param df Feature matrix with one row per sample and one column per mutation @@ -197,37 +584,55 @@ DLBCLone_train_test_plot = function(test_df, #' @param umap_out Optional UMAP output from a previous run. If provided, the function #' will use this model to project the data instead of re-running UMAP. This is useful #' for reproducibility and for using the same UMAP model on different datasets. -#' @param join_column The column name in the metadata data frame that contains the sample IDs (default sample_id). +#' @param core_features A vector of column names in df that should be multiplied by +#' core_feature_multiplier. These features are considered "core" and will be weighted +#' more heavily in the UMAP embedding. +#' @param core_feature_multiplier A numeric value that will be used to multiply the +#' core features specified in core_features. Default is 1.5. +#' @param hidden_features A vector of column names in df that should be dropped from the +#' UMAP input. These features are considered "hidden" and will not be used in the UMAP embedding. #' @param n_neighbors Passed to UMAP2. The number of neighbors to consider when calculating the UMAP embedding. #' @param min_dist Passed to UMAP2. The minimum distance between points in the UMAP embedding. #' @param metric Passed to UMAP2. The distance metric to use for calculating distances between points. #' @param n_epochs Passed to UMAP2. The number of epochs to run the UMAP algorithm. #' @param init Passed to UMAP2. The initialization method for the UMAP algorithm. +#' @param ret_model additional argument #' @param na_vals How to deal with NA values. Two options are "drop", which #' will remove all columns containing at least one NA or "to_zero", which sets #' all NA to zero and leaves the column intact. +#' @param join_column The column name in the metadata data frame that contains the sample IDs (default sample_id). #' @param seed Passed to UMAP2. The random seed for reproducibility. -#' @param ret_model additional argument +#' @param target_column The column in the metadata that contains the target variable for classification. +#' @param target_metric The distance metric to use for calculating distances to the target variable. +#' @param target_weight The weight to assign to the target variable in the UMAP embedding. Default is 0.5. +#' @param calc_dispersion If TRUE, calculates the dispersion of the UMAP embedding. +#' @param algorithm The UMAP algorithm to use. Default is "tumap", which uses the TUMAP implementation. +#' @param make_plot If TRUE, creates a plot of the UMAP embedding. #' -#' @import uwot +#' @import uwot dplyr tidyr tibble readr #' @export #' #' @examples +#' +#' library(GAMBLR.predict) #' -#' umap_outs = make_and_annotate_umap(df=gambl_train_lymphgen, -#' min_dist = 0, -#' n_neighbors = 55, -#' init="spectral", -#' n_epochs = 1500, -#' #seed=best_params$seed, -#' metadata=gambl_train_meta_dlbclass, -#' ret_model=T, -#' metric="cosine") -#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) #' make_and_annotate_umap = function(df, metadata, umap_out, + truth_column = "lymphgen", + core_features = NULL, + core_feature_multiplier = 1.5, + hidden_features = NULL, n_neighbors=55, min_dist=0, metric="cosine", @@ -236,108 +641,1481 @@ make_and_annotate_umap = function(df, ret_model=TRUE, na_vals = "drop", join_column="sample_id", - seed=42){ + seed=12345, + target_column, + target_metric="euclidean", + target_weight=0.5, + calc_dispersion = FALSE, + algorithm = "tumap", + make_plot = FALSE + ){ + + if(!missing(umap_out)){ + model_provided = TRUE + }else{ + model_provided = FALSE + } + if("sample_id" %in% colnames(df)){ + df = df %>% column_to_rownames(var = "sample_id") + } original_n = nrow(df) if(na_vals == "to_zero"){ df[is.na(df)] = 0 }else if(na_vals == "drop"){ - df <- df[, colSums(is.na(df)) == 0] + if(any(sapply(df, is.factor))){ + numeric_cols = names(df)[!sapply(df, is.factor)] + df <- df[, colSums(is.na(df[,numeric_cols])) == 0] + rs = rowSums(df[,numeric_cols],na.rm=TRUE) + dropped_rows = df[rs==0,] + df = df[rs>0,] + } else{ + df <- df[, colSums(is.na(df)) == 0] + rs = rowSums(df,na.rm=TRUE) + dropped_rows = df[rs==0,] + df = df[rs>0,] + } + } - rs = rowSums(df,na.rm=TRUE) - df = df[rs>0,] + + if(missing(df)){ stop("provide a data frame or matrix with one row for each sample and a numeric column for each mutation feature") } - if(missing(metadata)){ - stop("metadata is required and should contain a column sample_id that matches the row names of your mutation data frame") + if(!is.null(core_features)){ + #print(paste(core_features,collapse=",")) + if(!is.numeric(core_feature_multiplier)){ + stop("core_feature_multiplier must be a numeric value") + } + if(!all(core_features %in% colnames(df))){ + stop("core_features must be a vector of column names in df") + } + #multiply the core features by the multiplier + ncore = length(core_features) + message(paste0("multiplying ",ncore," core features by ",core_feature_multiplier)) + df[core_features] = df[core_features] * core_feature_multiplier } + if(!is.null(hidden_features)){ + if(!all(hidden_features %in% colnames(df))){ + stop("hidden_features must be a vector of column names in df") + } - df= df[rownames(df) %in% metadata[[join_column]],] + message(paste0("dropping ",length(hidden_features)," hidden features")) + df = df %>% select(-any_of(hidden_features)) + } + no_feat_samples = NULL + if(!missing(metadata)){ + if(nrow(df)< original_n){ + nrem = original_n-nrow(df) + pct_rem = round(nrem/original_n*100,2) + message(paste0("removed ",nrem," (",pct_rem,"%) rows from the data that had no features")) + no_feat_samples = rownames(dropped_rows) + + } + keep_rows = rownames(df)[rownames(df) %in% metadata[[join_column]]] + df= df[keep_rows,] + no_feat_metadata = filter(metadata,!!sym(join_column) %in% no_feat_samples) + metadata= filter(metadata,!!sym(join_column) %in% rownames(df)) + + message(paste("kept",nrow(metadata),"rows of the data that have features and match the metadata provided")) + } + if(missing(umap_out)){ - umap_out = umap2(df %>% as.matrix(), - n_neighbors = n_neighbors, - min_dist = min_dist, - metric = metric, - ret_model = ret_model, - n_epochs=n_epochs, - init=init, - seed = seed, - n_threads = 1) # possibly add rng_type = "deterministic" -#IMPORTANT: n_threads must not be changed because it will break reproducibility + if(missing(target_column)){ + if(algorithm == "umap"){ + message("running umap2") + umap_out = umap2(df %>% as.matrix(), + n_neighbors = n_neighbors, + min_dist = min_dist, + metric = metric, + ret_model = ret_model, + n_epochs=n_epochs, + a=1.8956, + b = 0.806, + approx_pow=TRUE, + init=init, + seed = seed, + n_threads = 1, + batch = TRUE, + n_sgd_threads = 1, + rng_type = "deterministic") # possibly add rng_type = "deterministic" + #IMPORTANT: n_threads must never be changed because it will break reproducibility + + }else if(algorithm == "tumap"){ + + message("running tumap") + X = df + umap_args = list(X=X, + n_neighbors = n_neighbors, + metric = metric, + ret_model = ret_model, + n_epochs=n_epochs, + init=init, + seed = seed, + n_threads = 1, + batch = TRUE, + n_sgd_threads = 1, + rng_type = "deterministic") + + umap_out = tumap(X, + n_neighbors = n_neighbors, + metric = metric, + ret_model = ret_model, + n_epochs=n_epochs, + init=init, + seed = seed, + n_threads = 1, + batch = TRUE, + n_sgd_threads = 1, + rng_type = "deterministic") + message("done!") + }else{ + stop("unsupported algorithm option") + } + }else{ + #supervised + if(missing(metadata)){ + stop("metadata must be provided for supervised UMAP") + } + metadata[[target_column]] = factor(metadata[[target_column]]) + + umap_out = umap2(df %>% as.matrix(), + n_neighbors = n_neighbors, + min_dist = min_dist, + metric = metric, + ret_model = ret_model, + n_epochs=n_epochs, + init=init, + seed = seed, + n_threads = 1, + y = metadata[[target_column]], + target_metric = target_metric, + target_weight = target_weight, + rng_type = "deterministic" + ) + #IMPORTANT: n_threads must never be changed because it will break reproducibility + } + }else{ - umap_out = umap_transform(X=df, - model=umap_out$model) - ret_model = FALSE + message("transforming each data point individually using the provided UMAP model. This will take some time.") + umap_df = data.frame() + for(sample in rownames(df)){ + this_row = df[sample,] + this_umap_df = umap_transform(X=this_row, + model=umap_out$model, + seed=seed, + batch = TRUE, + n_threads = 1, + n_sgd_threads = 1) + this_umap_df = as.data.frame(this_umap_df) + + umap_df = bind_rows(umap_df,this_umap_df) + } + message("done") + + } - if(ret_model){ - umap_df = as.data.frame(umap_out$embedding) %>% rownames_to_column(join_column) + + + if(model_provided){ + # model was generated here + #message("model given to function") + + umap_df = as.data.frame(umap_df) %>% rownames_to_column(var=join_column) + }else{ - umap_df = as.data.frame(umap_out) %>% rownames_to_column(join_column) + umap_df = as.data.frame(umap_out$embedding) %>% rownames_to_column(join_column) + } + if(!missing(metadata)){ + umap_df = left_join(umap_df,metadata,by=join_column) } - umap_df = left_join(umap_df,metadata) - results = list() + + + results[["df"]]=umap_df results[["features"]] = df + results[["dropped_rows"]] = no_feat_samples + if(!missing(metadata)){ + results[["total_samples_available"]] = nrow(metadata) + nrow(no_feat_metadata) + results[["sample_metadata_no_features"]] = no_feat_metadata + } + if(ret_model){ results[["model"]]= umap_out } + if(make_plot){ + ms = make_umap_scatterplot( + umap_df, + colour_by = truth_column, + title = "UMAP projection" + ) + print(ms) + results[["plot"]] = ms + } return(results) } +#' Process KNN Vote Strings and Scores for Classification +#' +#' This function processes the raw neighbor label strings and weighted vote scores from k-nearest neighbor (KNN) classification results. +#' It computes per-class neighbor counts, weighted scores, and identifies the top group by count and score for each sample. +#' The function also supports custom logic for handling the "Other" class, including vote multipliers and purity requirements. +#' +#' @param df Data frame containing kNN results, including columns with neighbor labels and weighted votes. +#' @param raw_col Name of the column containing the comma-separated neighbor labels (default: "label"). +#' @param group_labels Character vector of all possible class labels to consider (default: c("EZB", "MCD", "ST2", "BN2", "N1", "Other")). +#' @param vote_labels_col Name of the column containing the comma-separated neighbor labels for weighted votes (default: "vote_labels"). +#' @param k Number of neighbors used in kNN (required). +#' @param other_vote_multiplier Multiplier for the "Other" class when determining if a sample should be reclassified as "Other" (default: 2). +#' @param score_purity_requirement Minimum ratio of top group score to "Other" score to assign a sample to the top group (default: 1). +#' @param weighted_votes_col Name of the column containing the comma-separated weighted votes (default: "weighted_votes"). +#' +#' @return Data frame with additional columns for per-class neighbor counts, scores, top group assignments, and summary statistics for each sample. +#' +#' @details +#' - Computes the number of neighbors for each class and the sum of weighted votes per class. +#' - Identifies the top group by count and by weighted score, and applies custom logic for the "Other" class if present. +#' - Adds columns for counts, scores, top group, top group score, score ratios, and optimized group assignments. +#' - Designed for downstream use in DLBCLone and similar kNN-based classification workflows. +#' +#' @import dplyr tidyr stringr purrr +#' @export +#' +#' @examples +#' /dontrun{ +#' library(GAMBLR.predict) +#' +#' result <- process_votes(knn_output_df, k = 7) +#' } +#' +process_votes <- function(df, + raw_col = "label", + group_labels = c("EZB", "MCD", "ST2", "BN2", "N1", "Other"), + vote_labels_col = "vote_labels", + k, + other_vote_multiplier = 2, + score_purity_requirement = 1, + weighted_votes_col = "weighted_votes", + other_class = "Other") { # <--- NEW ARG + if(missing(k)){ + stop("k value is required") + } + score_thresh = 2 * k + + count_labels_in_string <- function(string, labels) { + tokens <- str_split(string, ",")[[1]] + map_int(labels, ~ sum(tokens == .x)) + } + + extract_weighted_scores <- function(label_str, vote_str, labels) { + lbls <- str_split(label_str, ",")[[1]] + votes <- as.numeric(str_split(vote_str, ",")[[1]]) + map_dbl(labels, ~ sum(votes[lbls == .x])) %>% + set_names(paste0(labels, "_score")) + } + + get_top_score_group <- function(label_str, vote_str, labels) { + lbls <- str_split(label_str, ",")[[1]] + votes <- as.numeric(str_split(vote_str, ",")[[1]]) + scores_by_label <- set_names(map_dbl(labels, ~ sum(votes[lbls == .x])), labels) + top <- names(scores_by_label)[which.max(scores_by_label)] + value <- scores_by_label[[top]] + list(top_score_group = top, top_group_score = value) + } + + df_out <- df %>% + mutate(.id = row_number()) %>% + rowwise() %>% + mutate( + counts = list( + set_names( + count_labels_in_string(.data[[raw_col]], group_labels), + paste0(group_labels, "_NN_count") + ) + ), + top_group = { + cnts <- count_labels_in_string(.data[[raw_col]], group_labels) + group_labels[which.max(cnts)] + }, + scores = list( + if (!is.null(vote_labels_col) && !is.null(weighted_votes_col)) { + extract_weighted_scores( + .data[[vote_labels_col]], + .data[[weighted_votes_col]], + group_labels + ) + } else { + set_names(rep(0, length(group_labels)), paste0(group_labels, "_score")) + } + ), + score_summary = list( + if (!is.null(vote_labels_col) && !is.null(weighted_votes_col)) { + get_top_score_group( + .data[[vote_labels_col]], + .data[[weighted_votes_col]], + group_labels + ) + } else { + list(top_score_group = NA_character_, top_group_score = NA_real_) + } + ) + ) %>% + ungroup() %>% + unnest_wider(counts) %>% + unnest_wider(scores) %>% + unnest_wider(score_summary) %>% + rowwise() %>% + mutate( + top_group_count = get(paste0(top_group, "_NN_count")) + ) %>% + ungroup() + + # Optional adjustments based on external columns (if present) + if (!is.null(other_class) && other_class %in% group_labels) { # <--- CONDITIONAL + if ("neighbors_other" %in% colnames(df)) { + df_out <- df_out %>% + mutate(!!sym(paste0(other_class, "_count")) := neighbors_other) + } + if ("other_weighted_votes" %in% colnames(df)) { + df_out <- df_out %>% + mutate(!!sym(paste0(other_class, "_score")) := other_weighted_votes) + } + if (all(c("top_group_count", paste0(other_class, "_count")) %in% colnames(df_out))) { + df_out <- df_out %>% + mutate(by_vote = top_group_count) %>% + mutate(by_vote_opt = ifelse(top_group_count * other_vote_multiplier > !!sym(paste0(other_class, "_count")), top_group, other_class)) + } + df_out <- mutate(df_out, + by_score = top_score_group, + score_ratio = top_group_score / !!sym(paste0(other_class, "_score")), + by_score_opt = ifelse(score_ratio > score_purity_requirement | top_group_score > score_thresh, top_score_group, other_class)) + } else { + # Fallback: if no "other" class exists, keep by_score/by_vote as top_group + df_out <- mutate(df_out, + by_score = top_score_group, + score_ratio = NA_real_, + by_score_opt = top_score_group, + by_vote_opt = top_group) + } + + return(df_out) +} + +#' Optimize Purity Threshold for Classification Assignment +#' +#' This function searches for the optimal purity threshold to assign samples to their predicted class or to "Other" based on the score ratio in processed kNN vote results. +#' It iteratively tests a range of purity thresholds, updating the predicted class if the score ratio meets or exceeds the threshold, and computes the accuracy for each threshold. +#' The function returns the best accuracy achieved and the corresponding purity threshold. +#' +#' @param optimized_model_object An object containing the optimized model parameters and predictions. +#' @param vote_df Data frame containing the kNN vote results, including columns for sample IDs, predicted labels, and scores. +#' @param mode The mode of operation, which determines the output column names (e.g., "DLBCLone_w"). +#' @param optimize_by The metric to optimize by, either "balanced_accuracy" +#' @param truth_column Name of the column in `processed_votes` containing the true class labels. +#' @param all_classes Vector of all possible class labels to consider (default: c("MCD", "EZB", "BN2", "N1", "ST2", "Other")). +#' @param k The number of neighbors used in the kNN classification. +#' @param exclude_other_for_accuracy Logical indicating whether to exclude the "Other" class from accuracy calculations (default: FALSE). +#' +#' @return A list with two elements: `best_accuracy` (numeric, the highest accuracy achieved) and `best_purity_threshold` (numeric, the threshold at which this accuracy was achieved). +#' +#' @details +#' - For each threshold in the range 0.1 to 0.95 (step 0.05), the function updates the prediction column to assign the class from `by_score_opt` if the score ratio meets the threshold, otherwise assigns "Other". +#' - Accuracy is computed as the proportion of correct assignments (diagonal of the confusion matrix). +#' - The function is intended for use in optimizing classification purity in kNN-based workflows, especially when distinguishing between confident class assignments and ambiguous ("Other") cases. +#' +#' @import dplyr tidyr caret +#' @export +#' +#' @examples +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' result <- optimize_purity(processed_votes, prediction_column = "pred_label", truth_column = "true_label") +#' } +#' +optimize_purity <- function(optimized_model_object, + vote_df, + mode, + optimize_by = "balanced_accuracy", #allowed: harmonic_mean, overall_accuracy + truth_column, + all_classes = c("MCD","EZB","BN2","N1","ST2","Other"), + k, + cap_classification_rate = 1, + exclude_other_for_accuracy = FALSE, + other_class = "Other") { # <--- NEW ARG + + out_column = "DLBCLone_wo" + + if(!missing(optimized_model_object)){ + if(!is.null(optimized_model_object$best_params)){ + if(!missing(k)){ + message("k is provided in the optimized_model_object, ignoring the k parameter") + } + k = optimized_model_object$best_params$k + }else{ + stop("optimized_model_object must contain best_params with k value") + } + vote_df = optimized_model_object$predictions + } + + some_classes <- if (!is.null(other_class) && other_class %in% all_classes) { + all_classes[all_classes != other_class] + } else { + all_classes + } + score_thresh = 2 * k + + processed_votes <- process_votes(vote_df, + group_labels = all_classes, + k = k, + score_purity_requirement = 0.5, + other_class = other_class) # <--- pass it through + + if(!truth_column %in% colnames(processed_votes)){ + stop("truth_column must be a column in processed_votes") + } + best_accuracy <- 0 + best_purity_threshold <- 0 + + processed_votes <- mutate(processed_votes, !!sym(truth_column) := as.character(!!sym(truth_column))) + + if (!is.null(other_class) && other_class %in% all_classes) { + no_other_df <- processed_votes %>% + filter(!!sym(truth_column) != other_class) + } else { + no_other_df <- processed_votes + } + + for(purity_threshold in seq(3, 0, -0.05)){ + updated_votes <- processed_votes %>% + mutate(!!sym(out_column) := if (!is.null(other_class) && other_class %in% all_classes) { + ifelse(score_ratio >= purity_threshold | top_group_score > score_thresh, by_score, other_class) + } else { + by_score + }) + + updated_no_other_df <- no_other_df %>% + mutate(!!sym(out_column) := if (!is.null(other_class) && other_class %in% all_classes) { + ifelse(score_ratio >= purity_threshold | top_group_score > score_thresh, by_score, other_class) + } else { + by_score + }) + + if(!exclude_other_for_accuracy || is.null(other_class) || !(other_class %in% all_classes)){ + xx <- select(updated_votes, sample_id, !!sym(truth_column), !!sym(out_column)) %>% + mutate(match = !!sym(truth_column) == !!sym(out_column)) %>% + group_by(match) %>% + summarise(concordant = sum(match == TRUE), discordant = sum(match == FALSE), .groups = "drop") %>% + summarise(all_conc = sum(concordant), dis = sum(discordant), total = all_conc + dis, percent = 100 * sum(concordant) / total) + + updated_votes <- mutate(updated_votes, !!sym(out_column) := factor(!!sym(out_column), levels = all_classes)) + updated_votes <- mutate(updated_votes, !!sym(truth_column) := factor(!!sym(truth_column), levels = all_classes)) + confusion_matrix <- table(updated_votes[[truth_column]], updated_votes[[out_column]]) + conf_matrix <- confusionMatrix(updated_votes[[out_column]], updated_votes[[truth_column]]) + } else { + xx <- select(updated_no_other_df, sample_id, !!sym(truth_column), !!sym(out_column)) %>% + filter(!!sym(truth_column) != other_class) %>% + mutate(match = !!sym(truth_column) == !!sym(out_column)) %>% + group_by(match) %>% + summarise(concordant = sum(match == TRUE), discordant = sum(match == FALSE), .groups = "drop") %>% + summarise(all_conc = sum(concordant), dis = sum(discordant), total = all_conc + dis, percent = 100 * sum(concordant) / total) + + updated_no_other_df <- mutate(updated_no_other_df, + !!sym(out_column) := factor(!!sym(out_column), levels = all_classes), + !!sym(truth_column) := factor(!!sym(truth_column), levels = all_classes)) + confusion_matrix <- table(updated_no_other_df[[truth_column]], updated_no_other_df[[out_column]]) + conf_matrix <- confusionMatrix(updated_no_other_df[[out_column]], updated_no_other_df[[truth_column]]) + } + + + acc_details <- report_accuracy(updated_votes, + truth = truth_column, + pred = out_column) + + bal_acc <- mean(conf_matrix$byClass[, "Balanced Accuracy"], na.rm = TRUE) + classification_rate = acc_details$classification_rate + if(classification_rate <= cap_classification_rate){ + if(optimize_by == "harmonic_mean"){ + accuracy = acc_details$harmonic_mean + }else if( optimize_by == "overall_accuracy"){ + accuracy = conf_matrix$overall[["Accuracy"]] + + }else{ + accuracy = bal_acc + } + + + if(accuracy > best_accuracy){ + best_accuracy <- accuracy + best_purity_threshold <- purity_threshold + } + }else{ + message(paste0("skipping purity threshold ",purity_threshold," because classification rate ",round(classification_rate,3), + " exceeds cap of ",cap_classification_rate)) + } + } + + best_out <- processed_votes %>% + mutate(!!sym(out_column) := if (!is.null(other_class) && other_class %in% all_classes) { + ifelse(score_ratio >= best_purity_threshold | top_group_score > score_thresh, by_score, other_class) + } else { + by_score + }) + + if(missing(optimized_model_object)){ + optimized_model_object <- list() + optimized_model_object$best_params <- list( + k = k, + purity_threshold = best_purity_threshold, + accuracy = best_accuracy, + num_classes = length(unique(best_out$predicted_label)), + num_features = ncol(best_out) - 3, + seed = 12345 + ) + } + optimized_model_object$predictions <- best_out + optimized_model_object$best_accuracy <- best_accuracy + optimized_model_object$best_purity_threshold <- best_purity_threshold + optimized_model_object$maximize <- optimize_by + optimized_model_object$score_thresh <- score_thresh + optimized_model_object$exclude_other_for_accuracy <- exclude_other_for_accuracy + + return(optimized_model_object) +} + +#' Predict DLBCLone Classes for New Samples Using a Trained KNN Model +#' +#' Applies a previously optimized DLBCLone KNN model to predict class labels for new (test) samples. +#' This function combines the training and test feature matrices, ensures feature compatibility, and uses the +#' parameters from a DLBCLone KNN optimization run to classify the test samples. Optionally, runs in iterative mode +#' for more stable results when predicting multiple samples. +#' +#' @param train_df Data frame or matrix of features for training samples (rows = samples, columns = features). +#' @param test_df Data frame or matrix of features for test samples to be classified. +#' @param metadata Data frame with metadata for all samples, including at least a \code{sample_id} column. +#' @param core_features Optional character vector of feature names to upweight in the KNN calculation. +#' @param core_feature_multiplier Numeric. Multiplier to apply to core features (default: 1.5). +#' @param hidden_features Optional character vector of feature names to exclude from the analysis. +#' @param DLBCLone_KNN_out List. Output from a previous call to \code{DLBCLone_KNN} containing optimized parameters. (Required) +#' @param mode Character. If \code{"iterative"}, runs KNN prediction for each test sample individually (recommended for stability). +#' +#' @return A list containing the KNN prediction results for the test samples, including predicted class labels and scores. +#' +#' @details +#' - Ensures that the feature columns in \code{train_df} and \code{test_df} are compatible. +#' - If \code{mode = "iterative"}, runs KNN prediction for each test sample one at a time. +#' - Uses the parameters (e.g., k, feature weights) from the provided \code{DLBCLone_KNN_out} object. +#' - Returns the same structure as \code{DLBCLone_KNN}, with predictions for the test samples. +#' +#' @import dplyr tidyr tibble readr +#' @export +#' +#' @examples +#' # Assuming you have run DLBCLone_KNN to get optimized parameters: +#' # model_out <- DLBCLone_KNN(train_features, train_metadata, ...) +#' # Predict on new samples: +#' +#' library(GAMBLR.predict) +#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' dlbcl_knn <- DLBCLone_KNN( +#' features_df = all_full_status, +#' metadata = dlbcl_meta +#' ) +#' +#' predictions <- DLBCLone_KNN_predict( +#' train_df = all_full_status, +#' test_df = lyseq_validation_data, +#' metadata = dlbcl_meta, +#' DLBCLone_KNN_out = dlbcl_knn, +#' mode = "batch" +#' ) +#' +#' +DLBCLone_KNN_predict <- function(train_df, + test_df, + metadata, + DLBCLone_KNN_out, + mode = "batch", + truth_column = "lymphgen", + other_class = "Other") { # <--- NEW ARG + if(missing(DLBCLone_KNN_out)){ + stop("DLBCLone_KNN_out must be provided, run DLBCLone_KNN first to get the optimal parameters") + } + nsamp = nrow(test_df) + message(paste0("Running DLBCLone KNN individually on ", nsamp, " samples")) + if(nsamp > 1){ + if(mode != "iterative"){ + warning("Running DLBCLone KNN on multiple samples at once is not recommended as the result may be unstable", + " and may not reflect the true classification of each sample. ", + "Running in iterative mode is recommended for more stable results.") + } + } + + if(any(!colnames(test_df) %in% colnames(train_df))){ + stop("test_df should not contain any features that are not in train_df. ", + "Please check the column names of test_df and train_df.") + } + combined_df = bind_rows(train_df, test_df) + if(any(!colnames(train_df) %in% colnames(test_df))){ + message("filling in missing features in test_df with zeros") + combined_df[is.na(combined_df)] = 0 + } + if(mode == "iterative"){ + predictions_list = list() + for(i in seq_len(nrow(test_df))){ + message("iteration:", i, "of", nrow(test_df)) + combined_df = bind_rows(train_df, test_df[i,]) + model_out = DLBCLone_KNN( + features_df = combined_df, + metadata = metadata, + DLBCLone_KNN_out = DLBCLone_KNN_out, + predict_unlabeled = TRUE, + other_class = other_class # <--- pass it in + ) + predictions_list[[i]] = model_out$unlabeled_predictions + } + all_predictions = do.call("bind_rows", predictions_list) + return(all_predictions) + } else { + model_out = DLBCLone_KNN( + features_df = combined_df, + metadata = metadata, + DLBCLone_KNN_out = DLBCLone_KNN_out, + truth_column = DLBCLone_KNN_out$truth_column, + truth_classes = DLBCLone_KNN_out$truth_classes, + predict_unlabeled = TRUE, + other_class = DLBCLone_KNN_out$other_class + ) + return(model_out) + } +} + +#' Run DLBCLone KNN Classification +#' +#' Weighted KNN on a feature (mutation) matrix with optional upweighting of +#' user-specified "core" features, optional exclusion of "hidden" features, +#' and optional optimization of an explicit outgroup (e.g. "Other"). +#' +#' This version removes hard-coded LymphGen class names and instead derives the +#' in-group classes and the outgroup column name from the arguments +#' \code{truth_classes} and \code{other_class}. It keeps backward compatibility +#' for the default LymphGen-like usage. +#' +#' @param features_df Numeric matrix/data.frame (rows = samples, cols = features). +#' Row names must be sample IDs. +#' @param metadata Data frame with at least \code{sample_id} and the ground-truth +#' label column given in \code{truth_column}. +#' @param core_features Character vector of feature names to upweight (optional). +#' @param core_feature_multiplier Numeric multiplier for \code{core_features}. +#' @param hidden_features Character vector of feature names to drop (optional). +#' @param min_k,max_k Integer K range to explore when optimizing. +#' @param truth_column Name of metadata column with ground-truth class labels. +#' @param truth_classes Character vector of all classes to consider (including +#' \code{other_class} if you intend to optimize for it). +#' @param other_class Name of the explicit outgroup class (default: "Other"). +#' @param optimize_for_other Logical; if TRUE, computes a separate "other" +#' score (ratio) and searches a purity threshold; if FALSE, treats all +#' classes symmetrically. +#' @param predict_unlabeled If TRUE, re-runs KNN to classify samples that were +#' present in \code{features_df} but not in \code{metadata}. +#' @param plot_samples Optional vector of sample_ids to keep in example plots. +#' @param DLBCLone_KNN_out Optional prior result; if supplied, its learned +#' parameters are reused (skip optimization). +#' @param seed Random seed. +#' @param epsilon Small value added to distances before weighting. +#' @param weighted_votes If FALSE, neighbors are unweighted (equal votes). +#' @param skip_umap If TRUE, skip layout optimization plots at the end. +#' +#' @return A list with fields including: +#' \item{predictions}{Per-sample vote/score summary and predicted labels} +#' \item{DLBCLone_k_best_k}{Best K found} +#' \item{DLBCLone_k_purity_threshold}{Best purity threshold (if applicable)} +#' \item{DLBCLone_k_accuracy}{Best accuracy metric achieved} +#' \item{truth_classes, truth_column}{Echoed arguments} +#' \item{unlabeled_predictions}{Predictions for unlabeled samples (if requested)} +#' \item{df}{Annotated layout for plotting (when built in this run)} +#' \item{plot_truth, plot_predicted}{ggplots when built in this run} +#' +#' @import dplyr tidyr uwot ggplot2 purrr caret +#' @export +#' +#' @examples +#' +#' library(GAMBLR.predict) +#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' dlbcl_knn <- DLBCLone_KNN( +#' features_df = all_full_status, +#' metadata = dlbcl_meta +#' ) +#' +DLBCLone_KNN <- function(features_df, + metadata, + core_features = NULL, + core_feature_multiplier = 1.5, + hidden_features = NULL, + min_k = 5, + max_k = 60, + truth_column = "lymphgen", + truth_classes = c("EZB", "BN2", "ST2", "MCD", "N1", "Other"), + other_class = "Other", + optimize_for_other = TRUE, + predict_unlabeled = FALSE, + plot_samples = NULL, + DLBCLone_KNN_out = NULL, + seed = 12345, + epsilon = 0.001, + weighted_votes = TRUE, + skip_umap = FALSE) { + + # In-group classes (exclude the explicit outgroup label if present) + class_levels <- setdiff(truth_classes, other_class) + + if(!missing(DLBCLone_KNN_out)){ + core_features <- DLBCLone_KNN_out$core_features + core_feature_multiplier <- DLBCLone_KNN_out$core_feature_multiplier + hidden_features <- DLBCLone_KNN_out$hidden_features + optimize_for_other <- DLBCLone_KNN_out$optimize_for_other + } + + # Upweight core features, drop hidden features + if(!is.null(core_features)){ + if(!is.numeric(core_feature_multiplier)){ + stop("core_feature_multiplier must be a numeric value") + } + if(!all(core_features %in% colnames(features_df))){ + stop("core_features must be a vector of column names in features_df") + } + ncore <- length(core_features) + message(paste0("multiplying ", ncore, " core features by ", core_feature_multiplier)) + features_df[, core_features] <- features_df[, core_features] * core_feature_multiplier + } + if(!is.null(hidden_features)){ + if(!all(hidden_features %in% colnames(features_df))){ + stop("hidden_features must be a vector of column names in features_df") + } + message(paste0("dropping ", length(hidden_features), " hidden features")) + features_df <- features_df %>% dplyr::select(-dplyr::any_of(hidden_features)) + } + + # Partition rows by presence in metadata and by empty feature rows + exclude_df <- features_df[!row.names(features_df) %in% metadata$sample_id, , drop = FALSE] + features_df <- features_df[ rownames(features_df) %in% metadata$sample_id, , drop = FALSE] + + df_empty <- features_df[rowSums(features_df) == 0, , drop = FALSE] + sample_metadata_no_features <- dplyr::filter(metadata, sample_id %in% rownames(df_empty)) + features_df <- features_df[rowSums(features_df) > 0, , drop = FALSE] + + exclude_empty <- exclude_df[rowSums(exclude_df) == 0, , drop = FALSE] + exclude_df <- exclude_df[rowSums(exclude_df) > 0, , drop = FALSE] + + if (is.null(DLBCLone_KNN_out)) { + # ------------------------ + # Optimize over K (and purity threshold if requested) + # ------------------------ + overall_best_acc_k <- 0 + overall_best_acc <- 0 + overall_best_thresh <- 0 + best_pred <- NULL + + metadata <- metadata %>% dplyr::filter(sample_id %in% rownames(features_df)) + metadata_simple <- metadata %>% dplyr::select(sample_id, !!rlang::sym(truth_column)) + + message("Finding all nearest neighbors up to k=", max_k + 1, " using cosine distance") + nn_u <- uwot::umap(features_df, + n_neighbors = max_k + 1, + ret_nn = TRUE, + metric = "cosine", + seed = seed, + n_threads = 1, + batch = TRUE, + n_sgd_threads = 1, + rng_type = "deterministic") + + fkn_ids <- nn_u$nn$cosine$idx + fkn_dists <- nn_u$nn$cosine$dist + rownames(fkn_dists) <- rownames(features_df) + rownames(fkn_ids) <- rownames(features_df) + + fkn_dists <- as.data.frame(fkn_dists) %>% dplyr::select(-1) + if (weighted_votes) { + fkn_weighted <- round(1 / (fkn_dists + epsilon), 7) + } else { + fkn_weighted <- fkn_dists + fkn_weighted[] <- 1 + } + + fkn_ids <- as.data.frame(fkn_ids) %>% dplyr::select(-1) + truth_index <- metadata[[truth_column]] + names(truth_index) <- metadata$sample_id + + fkn_ids_named <- apply(fkn_ids, 2, function(x) rownames(fkn_ids)[x]) + rownames(fkn_ids_named) <- rownames(fkn_ids) + + fkn_ids_truth <- apply(fkn_ids_named, 2, function(x) truth_index[x]) + rownames(fkn_ids_truth) <- rownames(fkn_ids) + fkn_ids_truth <- as.data.frame(fkn_ids_truth) + + fkn_ids_long <- tibble::rownames_to_column(fkn_ids_truth, "sample_id") %>% + tidyr::pivot_longer(-sample_id, names_to = "column", values_to = "class") + fkn_weighted_long <- tibble::rownames_to_column(fkn_weighted, "sample_id") %>% + tidyr::pivot_longer(-sample_id, names_to = "column", values_to = "vote") + + fkn_votes <- dplyr::left_join(fkn_ids_long, fkn_weighted_long, + by = c("sample_id", "column")) + + for (k in seq(max_k, min_k, by = -5)) { + message(paste0("Running DLBCLone KNN with k=", k)) + + fkn_weighted <- fkn_votes %>% + dplyr::group_by(sample_id) %>% + dplyr::slice_head(n = k) %>% + dplyr::ungroup() %>% + dplyr::group_by(sample_id, class) %>% + dplyr::summarize(weighted_vote = sum(vote), .groups = "drop") + + if (optimize_for_other) { + # Treat outgroup separately; compute top group among in-groups only + fkn_weighted <- fkn_weighted %>% + tidyr::pivot_wider( + id_cols = "sample_id", + names_from = "class", + values_from = "weighted_vote", + values_fill = 0 + ) %>% + dplyr::select(sample_id, dplyr::all_of(union(class_levels, other_class))) %>% + dplyr::rowwise() %>% + dplyr::mutate( + idx = which.max(dplyr::c_across(dplyr::all_of(class_levels))), + top_class = class_levels[idx], + top_class_count = dplyr::c_across(dplyr::all_of(class_levels))[idx] + ) %>% + dplyr::select(-idx) %>% + dplyr::ungroup() %>% + dplyr::mutate( + top_class = ifelse(top_class_count == 0, other_class, top_class), + "{other_class}" := round(.data[[other_class]], 4), + top_class_count = round(top_class_count, 4) + ) %>% + dplyr::rename( + "{other_class}_score" := dplyr::all_of(other_class) + ) %>% + dplyr::mutate( + by_score = top_class, + top_group_score = top_class_count, + score_ratio = round(top_group_score / .data[[paste0(other_class, "_score")]], 4) + ) + + } else { + # All classes symmetric; still provide an {other_class}_score column + fkn_weighted <- fkn_weighted %>% + tidyr::pivot_wider( + id_cols = "sample_id", + names_from = "class", + values_from = "weighted_vote", + values_fill = 0 + ) %>% + dplyr::select(sample_id, dplyr::all_of(union(class_levels, other_class))) %>% + dplyr::rowwise() %>% + dplyr::mutate( + idx = which.max(dplyr::c_across(dplyr::all_of(class_levels))), + top_class = class_levels[idx], + top_class_count = dplyr::c_across(dplyr::all_of(class_levels))[idx] + ) %>% + dplyr::select(-idx) %>% + dplyr::ungroup() %>% + dplyr::mutate( + top_class = ifelse(top_class_count == 0, other_class, top_class), + "{other_class}_score" := dplyr::if_else( + !is.na(.data[[other_class]]), .data[[other_class]], 0 + ) + ) %>% + dplyr::mutate( + by_score = top_class, + top_group_score = top_class_count, + score_ratio = 10 + ) + } + + # Attach truth labels for reporting + fkn_weighted <- dplyr::left_join( + fkn_weighted, + dplyr::select(metadata, sample_id, !!rlang::sym(truth_column)), + by = "sample_id" + ) + + score_thresh <- 1.5 * k + + # Grid search purity threshold if outgroup optimization is on + if (optimize_for_other) { + best <- 0 + best_thresh <- 0 + for (purity_threshold in seq(3, 0, -0.05)) { + updated_votes <- fkn_weighted %>% + dplyr::mutate(min_score = score_thresh) %>% + dplyr::mutate( + by_score_opt = ifelse( + score_ratio >= purity_threshold | top_group_score > score_thresh, + by_score, other_class + ) + ) + + acc <- report_accuracy(updated_votes, pred = "by_score_opt", truth = truth_column) + if (acc$mean_balanced_accuracy > best) { + best <- acc$mean_balanced_accuracy + best_thresh <- purity_threshold + } + } + } else { + updated_votes <- dplyr::mutate(fkn_weighted, by_score_opt = by_score) + acc <- report_accuracy(updated_votes, pred = "by_score_opt", truth = truth_column) + best <- acc$mean_balanced_accuracy + best_thresh <- 0 + } + + classes_for_span <- setdiff(intersect(class_levels, names(fkn_weighted)), other_class) + if (best > overall_best_acc) { + overall_best_acc <- best + overall_best_thresh <- best_thresh + overall_best_acc_k <- k + + updated_votes <- fkn_weighted %>% + dplyr::mutate(min_score = score_thresh) %>% + dplyr::mutate(DLBCLone_k = by_score) %>% + dplyr::mutate(DLBCLone_ko = ifelse( + score_ratio >= best_thresh | top_group_score > score_thresh, + by_score, other_class + )) %>% + dplyr::rowwise() %>% + dplyr::mutate( + valid_classes = { + scores <- dplyr::c_across(dplyr::all_of(classes_for_span)) + keep_idx <- which(replace(scores > score_thresh, is.na(scores), FALSE)) + if (length(keep_idx) == 0) other_class else paste(classes_for_span[keep_idx], collapse = "/") + } + ) %>% + dplyr::ungroup() + + best_pred <- updated_votes + message(paste0( + "New best accuracy: ", round(best, 3), + " at k=", k, " with purity threshold: ", round(best_thresh, 2) + )) + } else { + message(paste( + "best accuracy did not improve at k=", k, + " with purity threshold:", round(best_thresh, 2), + " accuracy:", round(best, 3) + )) + } + } # end for K + + } else { + message("Using DLBCLone_KNN_out provided, skipping KNN run") + best_pred <- DLBCLone_KNN_out$predictions + optimized_layout <- DLBCLone_KNN_out$df + overall_best_acc_k <- DLBCLone_KNN_out$DLBCLone_k_best_k + overall_best_acc <- DLBCLone_KNN_out$DLBCLone_k_accuracy + overall_best_thresh <- DLBCLone_KNN_out$DLBCLone_k_purity_threshold + } + + # ------------------------ + # Predict for unlabeled samples (optional) + # ------------------------ + unlabeled_predictions <- data.frame() + samples_no_metadata <- c( + rownames(exclude_df)[!rownames(exclude_df) %in% metadata$sample_id], + rownames(exclude_empty)[!rownames(exclude_empty) %in% metadata$sample_id] + ) + + if (predict_unlabeled && length(samples_no_metadata) > 0) { + message("Re-running KNN to include unlabeled samples. Will use K value and thresholds from optimized model, if provided") + message("will use newly provided features rather than recycling!") + + if (!is.null(DLBCLone_KNN_out)) { + k <- DLBCLone_KNN_out$DLBCLone_k_best_k + overall_best_thresh <- DLBCLone_KNN_out$DLBCLone_k_purity_threshold + best_thresh <- overall_best_thresh + } + + # Build an augmented metadata where the unlabeled have NA in the truth column + placeholder_metadata <- data.frame(sample_id = samples_no_metadata, stringsAsFactors = FALSE) + placeholder_metadata[[truth_column]] <- NA + metadata_merge <- dplyr::bind_rows(metadata, placeholder_metadata) + + k_buffer <- min(length(samples_no_metadata), overall_best_acc_k) + generous_k <- overall_best_acc_k + k_buffer + message("Finding all nearest neighbors up to k=", generous_k + 1, " using cosine distance") + + if (nrow(exclude_df) > 0) { + features_df_merge <- dplyr::bind_rows(features_df, exclude_df) + + nn_u <- uwot::umap(features_df_merge, + n_neighbors = generous_k + 1, + ret_nn = TRUE, + metric = "cosine", + seed = seed, + n_threads = 1, + batch = TRUE, + n_sgd_threads = 1, + rng_type = "deterministic") + fkn_ids <- nn_u$nn$cosine$idx + fkn_dists <- nn_u$nn$cosine$dist + rownames(fkn_dists) <- rownames(features_df_merge) + rownames(fkn_ids) <- rownames(features_df_merge) + + fkn_dists <- as.data.frame(fkn_dists) %>% dplyr::select(-1) + if (weighted_votes) { + fkn_weighted <- round(1 / (fkn_dists + epsilon), 7) + } else { + fkn_weighted <- fkn_dists + fkn_weighted[] <- 1 + } + + fkn_ids <- as.data.frame(fkn_ids) %>% dplyr::select(-1) + + truth_index <- metadata_merge[[truth_column]] + names(truth_index) <- metadata_merge$sample_id + + fkn_ids_named <- apply(fkn_ids, 2, function(x) rownames(fkn_ids)[x]) + rownames(fkn_ids_named) <- rownames(fkn_ids) + + fkn_ids_truth <- apply(fkn_ids_named, 2, function(x) truth_index[x]) + rownames(fkn_ids_truth) <- rownames(fkn_ids) + fkn_ids_truth <- as.data.frame(fkn_ids_truth) + + fk_neighbors_long <- as.data.frame(fkn_ids_named) %>% + tibble::rownames_to_column("sample_id") %>% + tidyr::pivot_longer(-sample_id, names_to = "column", values_to = "neighbor") + + fkn_ids_long <- tibble::rownames_to_column(fkn_ids_truth, "sample_id") %>% + tidyr::pivot_longer(-sample_id, names_to = "column", values_to = "class") + + fkn_weighted_long <- tibble::rownames_to_column(fkn_weighted, "sample_id") %>% + tidyr::pivot_longer(-sample_id, names_to = "column", values_to = "vote") + + fkn_votes <- dplyr::left_join(fkn_ids_long, fkn_weighted_long, by = c("sample_id", "column")) %>% + dplyr::left_join(fk_neighbors_long, by = c("sample_id", "column")) %>% + dplyr::filter(!is.na(class)) + + k <- overall_best_acc_k + + fkn_weighted <- fkn_votes %>% + dplyr::group_by(sample_id) %>% + dplyr::slice_head(n = k) %>% + dplyr::mutate(neighbor_id = paste(neighbor, collapse = ",")) %>% + dplyr::ungroup() + + fkn_weighted_neighbors <- fkn_weighted %>% + dplyr::select(sample_id, neighbor_id) %>% + dplyr::distinct() %>% + dplyr::filter(sample_id %in% rownames(exclude_df)) + + fkn_weighted <- fkn_weighted %>% + dplyr::group_by(sample_id, class) %>% + dplyr::summarize(weighted_vote = sum(vote), .groups = "drop") + + if (optimize_for_other) { + fkn_weighted <- fkn_weighted %>% + tidyr::pivot_wider( + id_cols = "sample_id", + names_from = "class", + values_from = "weighted_vote", + values_fill = 0 + ) %>% + dplyr::select(sample_id, dplyr::all_of(union(class_levels, other_class))) %>% + dplyr::rowwise() %>% + dplyr::mutate( + idx = which.max(dplyr::c_across(dplyr::all_of(class_levels))), + top_class = class_levels[idx], + top_class_count = dplyr::c_across(dplyr::all_of(class_levels))[idx] + ) %>% + dplyr::select(-idx) %>% + dplyr::ungroup() %>% + dplyr::mutate( + top_class = ifelse(top_class_count == 0, other_class, top_class), + "{other_class}" := round(.data[[other_class]], 4), + top_class_count = round(top_class_count, 4) + ) %>% + dplyr::rename( + "{other_class}_score" := dplyr::all_of(other_class) + ) %>% + dplyr::mutate( + by_score = top_class, + top_group_score = top_class_count, + score_ratio = round(top_group_score / .data[[paste0(other_class, "_score")]], 4) + ) + } else { + fkn_weighted <- fkn_weighted %>% + tidyr::pivot_wider( + id_cols = "sample_id", + names_from = "class", + values_from = "weighted_vote", + values_fill = 0 + ) %>% + dplyr::select(sample_id, dplyr::all_of(union(class_levels, other_class))) %>% + dplyr::rowwise() %>% + dplyr::mutate( + idx = which.max(dplyr::c_across(dplyr::all_of(class_levels))), + top_class = class_levels[idx], + top_class_count = dplyr::c_across(dplyr::all_of(class_levels))[idx] + ) %>% + dplyr::select(-idx) %>% + dplyr::ungroup() %>% + dplyr::mutate( + top_class = ifelse(top_class_count == 0, other_class, top_class), + "{other_class}_score" := dplyr::if_else( + !is.na(.data[[other_class]]), .data[[other_class]], 0 + ) + ) %>% + dplyr::mutate( + by_score = top_class, + top_group_score = top_class_count, + score_ratio = 10 + ) + } + + fkn_weighted <- fkn_weighted %>% + dplyr::filter(sample_id %in% rownames(exclude_df)) + + fkn_weighted <- dplyr::left_join( + fkn_weighted, + dplyr::select(metadata_merge, sample_id, !!rlang::sym(truth_column)), + by = "sample_id" + ) + + score_thresh <- 1.5 * k + + if (optimize_for_other) { + unlabeled_predictions <- fkn_weighted %>% + dplyr::mutate(DLBCLone_k = by_score) %>% + dplyr::mutate(DLBCLone_ko = ifelse( + score_ratio >= overall_best_thresh | top_group_score > score_thresh, + by_score, other_class + )) %>% + dplyr::left_join(., fkn_weighted_neighbors, by = "sample_id") + } else { + classes_for_span <- intersect(class_levels, names(fkn_weighted)) + other_score_col <- paste0(other_class, "_score") + + unlabeled_predictions <- fkn_weighted %>% + dplyr::mutate(DLBCLone_k = by_score) %>% + dplyr::mutate(DLBCLone_ko = by_score) %>% + dplyr::rowwise() %>% + dplyr::mutate( + valid_classes = { + scores <- c_across(all_of(classes_for_span)) + keep_idx <- which(replace(scores > .data[[other_score_col]], is.na(scores), FALSE)) + if (length(keep_idx) == 0) other_class else paste(classes_for_span[keep_idx], collapse = "/") + } + ) %>% + dplyr::ungroup() %>% + dplyr::left_join(., fkn_weighted_neighbors, by = "sample_id") + } + } else { + other_score_col <- paste0(other_class, "_score") + unlabeled_predictions <- data.frame(sample_id = rownames(exclude_empty), stringsAsFactors = FALSE) + unlabeled_predictions[[truth_column]] <- NA + unlabeled_predictions <- unlabeled_predictions %>% + dplyr::mutate( + DLBCLone_k = other_class, + DLBCLone_ko = other_class, + by_score = NA, + top_group_score = NA, + score_ratio = NA + ) + # Ensure truth_classes columns exist as NA + for (cls in truth_classes) { + if (!cls %in% colnames(unlabeled_predictions)) { + unlabeled_predictions[[cls]] <- NA + } + } + unlabeled_predictions[[other_score_col]] <- NA + # no neighbor_id information in this branch + } + } + + # ------------------------ + # Assemble return object + # ------------------------ + if (is.null(DLBCLone_KNN_out)) { + + if (nrow(sample_metadata_no_features) > 0) { + sample_metadata_no_features <- sample_metadata_no_features %>% + dplyr::select(sample_id, !!rlang::sym(truth_column)) %>% + dplyr::mutate( + DLBCLone_k = NA, + DLBCLone_ko = other_class, + by_score = NA, + top_group_score = NA, + score_ratio = NA + ) + best_pred <- dplyr::bind_rows(best_pred, sample_metadata_no_features) + } + + format_for_output <- function(x) { + x_rounded <- round(x, 4) + ifelse( + is.na(x_rounded), + NA_character_, + ifelse( + x_rounded == 0, + "0", + sub("\\.?0+$", "", format(x_rounded, scientific = FALSE, trim = TRUE, nsmall = 0)) + ) + ) + } + + to_return <- list( + predictions = best_pred %>% + dplyr::mutate(dplyr::across(where(is.numeric), ~ format_for_output(.))), + DLBCLone_k_best_k = overall_best_acc_k, + DLBCLone_k_purity_threshold = overall_best_thresh, + DLBCLone_k_accuracy = overall_best_acc, + truth_classes = truth_classes, + truth_column = truth_column, + sample_metadata_no_features = sample_metadata_no_features, + core_feature_multiplier = core_feature_multiplier, + core_features = core_features, + hidden_features = hidden_features, + seed = seed, + optimize_for_other = optimize_for_other, + unlabeled_predictions = NULL + ) + } else { + to_return <- DLBCLone_KNN_out + } + + # Fill in neighbors / predictions for unlabeled (if requested) + to_return$unlabeled_neighbors <- NULL + if (predict_unlabeled && length(samples_no_metadata) > 0) { + if (exists("fkn_weighted_neighbors")) { + to_return$unlabeled_neighbors <- fkn_weighted_neighbors %>% + tidyr::separate(neighbor_id, into = paste0("N", c(1:k)), sep = ",") + } + to_return$unlabeled_predictions <- unlabeled_predictions + + if (is.null(DLBCLone_KNN_out) & !skip_umap) { + message("Optimizing graph layout for visualization, predict_unlabeled = TRUE") + df_show <- dplyr::bind_rows(features_df, exclude_df) + + optimized <- make_and_annotate_umap(df_show, metadata = metadata)$df + optimized <- optimized %>% + dplyr::left_join( + dplyr::bind_rows( + dplyr::select(best_pred, -!!rlang::sym(truth_column)), + dplyr::select(unlabeled_predictions, -!!rlang::sym(truth_column)) + ), + by = "sample_id" + ) %>% + dplyr::mutate(!!rlang::sym(truth_column) := ifelse(is.na(.data[[truth_column]]), DLBCLone_ko, .data[[truth_column]])) + } else { + optimized <- if (!is.null(DLBCLone_KNN_out)) DLBCLone_KNN_out$features_df else NULL + } + } else { + if (is.null(DLBCLone_KNN_out)) { + message("Optimizing graph layout for visualization") + optimized <- make_and_annotate_umap(features_df, metadata = metadata)$df + optimized <- optimized %>% + dplyr::left_join(dplyr::select(best_pred, -!!rlang::sym(truth_column)), by = "sample_id") %>% + dplyr::mutate(!!rlang::sym(truth_column) := ifelse(is.na(.data[[truth_column]]), DLBCLone_ko, .data[[truth_column]])) + } else { + message("Using DLBCLone_KNN_out provided, skipping graph layout optimization") + optimized <- DLBCLone_KNN_out$features_df + } + } + + if(!is.null(optimized) & is.null(DLBCLone_KNN_out)){ + to_return$plot_truth <- basic_umap_scatterplot(optimized, plot_samples, colour_by = truth_column) + to_return$plot_predicted <- basic_umap_scatterplot(optimized, plot_samples, colour_by = "DLBCLone_ko") + to_return$df <- optimized + } + + # Ensure features_df is included in the return for downstream plotting + if (predict_unlabeled && length(samples_no_metadata) > 0) { + to_return$features_df <- if (exists("features_df_merge")) features_df_merge else features_df + } else { + to_return$features_df <- features_df + } + + to_return$type <- "DLBCLone_KNN" + to_return$pred_column <- "DLBCLone_ko" + to_return$other_class = other_class + return(to_return) +} + #' Optimize parameters for classifying samples using UMAP and k-nearest neighbor #' -#' @param combined_mutation_status_df Data frame with one row per sample and one column per mutation -#' @param metadata_df Data frame of metadata with one row per sample and three required columns: sample_id, dataset and lymphgen -#' @param truth_classes Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2","Other") -#' @param eval_group Specify whether certain rows will be evaluated and held out from training rather than using all samples. -#' @param umap_out additional argument -#' @param min_k additional argument -#' @param max_k additional argument -#' @param verbose additional argument -#' @param seed additional argument +#' @param combined_mutation_status_df Data frame with one row per sample and +#' one column per mutation +#' @param metadata_df Data frame of metadata with one row per sample and +#' three required columns: sample_id, dataset and lymphgen +#' @param umap_out The output of a previous run of make_and_annotate_umap. +#' If provided, the function will use this model to project the data +#' instead of re-running UMAP. +#' @param truth_classes Vector of classes to use for training and testing. +#' Default: c("EZB","MCD","ST2","N1","BN2","Other") +#' @param truth_column The column in the metadata_df that contains the true class labels. +#' Default: "lymphgen". +#' @param optimize_for_other Set to TRUE to optimize the threshold for +#' classifying samples as "Other" based on the relative proportion of +#' samples near the sample in UMAP space with the "Other" label. Rather than +#' treating Other as just another class, this will optimize the threshold for +#' a separate score that considers how many Other and non-Other samples are +#' in the neighborhood of the sample in question. This parameter will NOT change +#' the value in predicted_label. Instead, the predicted_label_optimized column +#' will contain the optimized label. Default: FALSE +#' @param eval_group If desired, use this to specify which rows will be +#' evaluated and held out from training rather than using all samples. +#' NOTE: this parameter will probably become deprecated! +#' @param min_k Starting k for knn (Default: 3) +#' @param max_k Ending k for knn (Default: 33) +#' @param verbose Whether to print verbose outputs to console +#' @param seed Random seed to use for reproducibility (default: 12345) +#' @param maximize Metric to use for optimization. Either "sensitivity" +#' (average sensitivity across all classes), "accuracy" +#' (actual accuracy across all samples) or "balanced_accuracy" (the mean of the +#' balanced accuracy values across all classes). Default: "balanced_accuracy" +#' @param exclude_other_for_accuracy Set to TRUE to exclude the +#' "Other" class from the 'lymphgen' column when calculating accuracy metrics +#' (passed to DLBCLone_optimize_params). Default: FALSE +#' @param weights_opt Vector of TRUE/FALSE values to indicate whether to use +#' weighted knn or not. Default: c(TRUE). #' #' @returns List of data frames with the results of the parameter optimization -#' including the best model, the associated knn parameters and the annotated UMAP output +#' including the best model, the associated knn parameters and the annotated +#' UMAP output as a data frame. The list also includes the predictions for the +#' "Other" class if it was included in the training and testing. +#' +#' @import dplyr tidyr tibble #' @export #' #' @examples #' -#' lymphgen_A53_DLBCLone = DLBCLone_optimize_params( -#' lgen_feat_status, #our binary feature matrix -#' a53_meta, #our metadata -#' umap_out = lymphgen_A53_all_feat_gambl, # force use existing UMAP fit -#' eval_group = NULL, # use all samples for evaluating accuracy -#' truth_classes = c("MCD","EZB","BN2","ST2","N1","A53","Other")) - +#' library(GAMBLR.predict) +#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) +#' +#' # Aim to maximize classification of samples into non-Other class +#' optimize_params_F <- DLBCLone_optimize_params( +#' all_full_status, +#' dlbcl_meta, +#' umap_out = make_umap, +#' truth_classes = c("MCD","EZB","BN2","N1","ST2","Other"), +#' optimize_for_other = F, +#' min_k=5, +#' max_k=23 +#' ) +#' +#' # Aim to maximize balanced accuracy while allowing samples to be +#' # unclassified (assigned to "Other") +#' +#' optimize_params_T <- DLBCLone_optimize_params( +#' all_full_status, +#' dlbcl_meta, +#' umap_out = make_umap, +#' truth_classes = c("MCD","EZB","BN2","N1","ST2","Other"), +#' optimize_for_other = T, +#' min_k=5, +#' max_k=23 +#' ) +#' DLBCLone_optimize_params = function(combined_mutation_status_df, metadata_df, umap_out, - truth_classes = c("EZB","MCD","ST2","N1","BN2","Other"), - eval_group = "Lacy", + truth_classes = c("EZB", + "MCD", + "ST2", + "N1", + "BN2", + "Other"), + truth_column = "lymphgen", + optimize_for_other = FALSE, + + eval_group = NULL, min_k=3, - max_k=30, + max_k=23, verbose = FALSE, - seed = 12345){ + seed = 12345, + maximize = "balanced_accuracy", #or "harmonic_mean" or "accuracy" + cap_classification_rate = 0.9, + exclude_other_for_accuracy = FALSE, + weights_opt = c(TRUE) + ) { + macro_f1 <- function(truth, pred, drop_other = TRUE, other_label = "Other", na_rm = TRUE) { + stopifnot(length(truth) == length(pred)) + truth <- factor(truth) + pred <- factor(pred, levels = levels(truth)) # align levels + + classes <- levels(truth) + if (drop_other && other_label %in% classes) { + classes <- setdiff(classes, other_label) + } + + f1_per_class <- sapply(classes, function(cls) { + tp <- sum(truth == cls & pred == cls) + fp <- sum(truth != cls & pred == cls) + fn <- sum(truth == cls & pred != cls) + + prec <- if ((tp + fp) == 0) NA_real_ else tp / (tp + fp) + rec <- if ((tp + fn) == 0) NA_real_ else tp / (tp + fn) + + if (is.na(prec) || is.na(rec) || (prec + rec) == 0) { + if (na_rm) return(NA_real_) else return(0) # choose policy + } + 2 * prec * rec / (prec + rec) + }) + + mean(f1_per_class, na.rm = na_rm) + + } + + exclude_other_for_accuracy = TRUE + na_opt = c("drop") num_class = length(truth_classes) - weights_opt = c(TRUE,FALSE) - threshs = seq(0,0.95,0.05) - #ks = seq(5,25,1) - ks = seq(min_k,max_k,1) + + threshs = seq(0,0.9,0.1) + + + ks = seq(min_k,max_k,2) results <- data.frame() best_params <- data.frame() - #best so far: - threshold = 0.3 - k = 6 + use_w = TRUE best_acc = 0 + best_acc_w = 0 + best_pred_w = NULL + best_k_w = 0 best_fit = NULL - + this_accuracy = 0 best_pred = NULL other_pred = NULL + best_w_score_thresh = NULL + best_w_purity = NULL + w_best_pred = NULL + ignore_self = TRUE + for(na_option in na_opt){ if(missing(umap_out)){ outs = make_and_annotate_umap(df=combined_mutation_status_df, @@ -352,7 +2130,10 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, na_vals = na_option) }else{ #project onto existing model instead of re-running UMAP - outs = make_and_annotate_umap(df=combined_mutation_status_df, + if(!missing(combined_mutation_status_df)){ + message("ignoring mutation status data frame and using features from umap_out instead") + } + outs = make_and_annotate_umap(df=umap_out$features, umap_out = umap_out, min_dist = 0, n_neighbors = 55, @@ -364,59 +2145,115 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, join_column="sample_id", na_vals = na_option) } - ignore_top = FALSE - if(is.null(eval_group)){ - ignore_top = TRUE - } + + for(use_w in weights_opt){ for(k in ks){ - - for (threshold in threshs){ - if(is.null(eval_group)){ - test_coords = filter(outs$df,lymphgen %in% truth_classes) %>% select(V1,V2) - train_coords = filter(outs$df,lymphgen %in% truth_classes) %>% select(V1,V2) - train_labels = filter(outs$df,lymphgen %in% truth_classes) %>% pull(lymphgen) + best_w_acc = NULL + message(paste("K:",k)) + test_coords = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% select(V1,V2) + train_coords = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% select(V1,V2) + train_labels = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% pull(!!sym(truth_column)) + train_ids = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% pull(sample_id) + test_ids = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% pull(sample_id) + rownames(train_coords) = train_ids + rownames(test_coords) = test_ids - }else{ - test_coords = filter(outs$df,dataset == eval_group) %>% select(V1,V2) - train_coords = filter(outs$df,dataset != eval_group,lymphgen %in% truth_classes) %>% select(V1,V2) - train_labels = filter(outs$df,dataset != eval_group,lymphgen %in% truth_classes) %>% pull(lymphgen) - } + + pred = weighted_knn_predict_with_conf( + train_coords = train_coords, + train_labels = train_labels, + test_coords = test_coords, + k=k, + conf_threshold =threshold, + other_class="Other", + use_weights = use_w, + ignore_self = ignore_self, + verbose = verbose) + + for(threshold in threshs){ if(verbose){ print(paste("k:",k,"threshold:",threshold,"use_weights:",use_w,"na_option:",na_option)) } + pred_thresh = mutate(pred, + predicted_label = ifelse(confidence >= threshold, predicted_label, "Other")) - pred = weighted_knn_predict_with_conf( - train_coords = train_coords, - train_labels = train_labels, - test_coords = test_coords, - k=k, - conf_threshold =threshold, - na_label="Other", - use_weights = use_w, - ignore_top = ignore_top, - verbose = verbose) - - if(is.null(eval_group)){ - xx_d = bind_cols(filter(outs$df,lymphgen %in% truth_classes) ,pred) - train_d = filter(outs$df,lymphgen %in% truth_classes) - }else{ - xx_d = bind_cols(filter(outs$df,dataset == eval_group) ,pred) - train_d = filter(outs$df,dataset != eval_group,lymphgen %in% truth_classes) + if(is.null(best_w_acc)){ + #DLBCLone_wo + some_outs = filter(outs$df,!!sym(truth_column) %in% truth_classes) + pred_with_truth = bind_cols(some_outs ,pred) + if(verbose){ + print(dim(pred_with_truth)) + print(table(pred_with_truth[[truth_column]])) + } + + pred_w = + optimize_purity( + vote_df = pred_with_truth, + mode = "DLBCLone_w", + truth_column = truth_column, + all_classes = truth_classes, + optimize_by = maximize, + cap_classification_rate = cap_classification_rate, + k = k + ) + + if(verbose){ + print("running optimize_purity for the first time at k:", k) + } + reported_accuracy = report_accuracy(pred_w$predictions,pred = "DLBCLone_wo") + best_w_acc = pred_w$best_accuracy + if(best_w_acc > best_acc_w){ + + best_acc_w = best_w_acc + best_k_w = k + best_fit = pred_with_truth + best_pred_w = pred_w$predictions + best_w_purity = pred_w$best_purity_threshold + best_w_score_thresh = pred_w$score_thresh + #print(head(pred_w$predictions)) + + conc = reported_accuracy$accuracy_no_other + f1 = macro_f1(pred_w$predictions[[truth_column]], + pred_w$predictions$DLBCLone_wo, + drop_other = T) + new_f1 = reported_accuracy$macro_f1 + harmonic_mean = reported_accuracy$harmonic_mean + mba = reported_accuracy$mean_balanced_accuracy + cr = reported_accuracy$classification_rate + + message("new best DLBCLone_wo:") + message(paste( + "Concordance:", + round(conc,3), + "F1:", + round(f1,5), + "MBA:", + round(best_w_acc,3), + "Classification rate:", + round(cr,5), + "Harmonic mean:", + round(harmonic_mean,5), + "purity:", + pred_w$best_purity_threshold, + "Score:", + best_w_score_thresh + )) + } + best_pred_w = pred_w$predictions + } + + train_d = filter(outs$df,!!sym(truth_column) %in% truth_classes) + xx_d = bind_cols(train_d ,pred_thresh) + if("Other" %in% truth_classes){ - xx_d$lymphgen = factor(xx_d$lymphgen) + xx_d[[truth_column]] = factor(xx_d[[truth_column]]) }else{ - if(is.null(eval_group)){ - test_coords = filter(outs$df,lymphgen %in% "Other") %>% select(V1,V2) - train_coords = filter(outs$df,lymphgen %in% truth_classes) %>% select(V1,V2) - train_labels = filter(outs$df,lymphgen %in% truth_classes) %>% pull(lymphgen) - }else{ - test_coords = filter(outs$df,dataset == eval_group,lymphgen %in% "Other") %>% select(V1,V2) - train_coords = filter(outs$df,dataset == eval_group,lymphgen %in% unique(c("Other",truth_classes))) %>% select(V1,V2) - train_labels = filter(outs$df,dataset == eval_group,lymphgen %in% unique(c("Other",truth_classes))) %>% pull(lymphgen) - } + test_coords = filter(outs$df,!!sym(truth_column) %in% "Other" | is.na(!!sym(truth_column))) %>% select(V1,V2) + train_coords = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% select(V1,V2) + train_labels = filter(outs$df,!!sym(truth_column) %in% truth_classes) %>% pull(!!sym(truth_column)) n_other = nrow(test_coords) if(!"Other" %in% truth_classes && n_other > 0){ @@ -426,22 +2263,20 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, test_coords = test_coords, k=k, conf_threshold =threshold, - na_label="Other", + other_class="Other", use_weights = use_w, - ignore_top = ignore_top, + ignore_self = ignore_self, verbose = verbose) - if(!is.null(eval_group)){ - xx_o = bind_cols(filter(outs$df,dataset == eval_group,lymphgen =="Other") ,pred_other) - }else{ - xx_o = bind_cols(filter(outs$df,lymphgen == "Other") ,pred_other) - } + xx_o = bind_cols(filter(outs$df,!!sym(truth_column) == "Other" | is.na(!!sym(truth_column))) ,pred_other) + } - xx_d$lymphgen = factor(xx_d$lymphgen,levels = c(unique(xx_d$lymphgen),"Other")) + xx_d[[truth_column]] = factor(xx_d[[truth_column]],levels = c(unique(xx_d[[truth_column]]),"Other")) } + true_factor = xx_d[[truth_column]] + pred_factor = factor(xx_d$predicted_label,levels = levels(xx_d[[truth_column]])) - true_factor = factor(xx_d$lymphgen,levels = levels(xx_d$lymphgen)) if(verbose){ print("true_factor") print(levels(true_factor )) @@ -449,23 +2284,71 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, print(levels(xx_d$predicted_label)) } - conf_matrix <- confusionMatrix(xx_d$predicted_label, true_factor) - + conf_matrix <- confusionMatrix(pred_factor, true_factor) + bal_acc <- conf_matrix$byClass[, "Balanced Accuracy"] # one per class sn <- conf_matrix$byClass[, "Sensitivity"] # one per class + if(verbose){ print(bal_acc) + print(conf_matrix$overall) + print(sn) } + + overall_accuracy <- conf_matrix$overall[["Accuracy"]] + + overall_sensitivity<- mean(sn[!names(sn) == "Class: Other"], na.rm = TRUE) + + if(optimize_for_other){ + + optimized_accuracy_and_thresh = optimize_outgroup(pred_factor, + true_factor, + xx_d$other_score, + all_classes = truth_classes, + maximize = maximize, + cap_classification_rate = cap_classification_rate, + exclude_other_for_accuracy = exclude_other_for_accuracy) + + out_opt_thresh = optimized_accuracy_and_thresh$threshold + if(maximize=="harmonic_mean"){ + out_opt_acc = optimized_accuracy_and_thresh$harmonic_mean + }else{ + out_opt_acc = optimized_accuracy_and_thresh$average_accuracy + } + - overall_balanced_accuracy <- mean(bal_acc, na.rm = TRUE) + }else{ + out_opt_acc = 0 + out_opt_thresh = 0 + } + if(exclude_other_for_accuracy){ + mean_balanced_accuracy = mean(bal_acc[!names(bal_acc) == "Class: Other"], na.rm = TRUE) + }else{ + mean_balanced_accuracy = mean(bal_acc) + } + if(maximize == "sensitivity"){ + this_accuracy = overall_sensitivity + }else if(maximize == "balanced_accuracy"){ + + this_accuracy =mean_balanced_accuracy + }else{ + + this_accuracy = overall_accuracy + } + if(out_opt_acc > this_accuracy){ + + this_accuracy = out_opt_acc + } row <- data.frame(k = k, threshold = threshold, use_weights = use_w, - accuracy = overall_balanced_accuracy, + optimized_accuracy = this_accuracy, + overall_accuracy = overall_accuracy, + mean_balanced_accuracy = mean_balanced_accuracy, + sensitivity = overall_sensitivity, N1_sn = unname(sn["Class: N1"]), BN2_sn= unname(sn["Class: BN2"]), ST2_sn = unname(sn["Class: ST2"]), - N1_bacc = unname(bal_acc["Class: N1"]), BN2_bacc = unname(bal_acc["Class: BN2"]), MCD_bacc = unname(bal_acc["Class: MCD"]), EZB_bacc = unname(bal_acc["Class: EZB"]), @@ -476,10 +2359,25 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, C3_bacc = unname(bal_acc["Class: C3"]), C4_bacc = unname(bal_acc["Class: C4"]), C5_bacc = unname(bal_acc["Class: C5"]), - ST2_bacc = unname(bal_acc["Class: ST2"]), na_option= na_option) - if(overall_balanced_accuracy > best_acc){ - best_acc = overall_balanced_accuracy - print(paste("best accuracy:",best_acc,"k:",k,"threshold:",threshold,"na:",na_option,"Balanced accuracy:",overall_balanced_accuracy)) + ST2_bacc = unname(bal_acc["Class: ST2"]), + threshold_outgroup = out_opt_thresh, + accuracy_out = out_opt_acc, + na_option= na_option) + + if(this_accuracy > best_acc){ + best_acc = this_accuracy + + xx_d = mutate(xx_d, predicted_label_optimized = ifelse(other_score > out_opt_thresh, "Other", predicted_label)) + no_other_acc = report_accuracy(xx_d,pred = "predicted_label_optimized")$mean_balanced_accuracy + message("new best accuracy DLBCLone_io:") + message(paste( + best_acc, + "without other: ", + no_other_acc, + "sensitivity:", + overall_sensitivity, + "threshold_outgroup:", + row$threshold_outgroup)) best_fit = outs best_pred = xx_d @@ -495,38 +2393,101 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, } } best_params$num_classes = num_class - best_params$num_features = ncol(best_fit$features) + best_params$num_features = ncol(umap_out$features) best_params$seed = seed test_coords = outs$df %>% select(V1,V2) - train_coords = filter(outs$df,lymphgen %in% truth_classes) %>% select(V1,V2) - train_labels = filter(outs$df,lymphgen %in% truth_classes) %>% pull(lymphgen) + + train_coords = outs$df %>% select(V1,V2) + train_ids = outs$df %>% pull(sample_id) + rownames(train_coords) = train_ids + train_labels = outs$df %>% pull(!!sym(truth_column)) + rownames(test_coords) = outs$df %>% pull(sample_id) + + if(verbose){ + print(paste("TOP score threshold:",best_w_score_thresh, "purity:", best_w_purity)) + } + pred = weighted_knn_predict_with_conf( train_coords = train_coords, train_labels = train_labels, test_coords = test_coords, - k=best_params$k, - conf_threshold =best_params$threshold, - na_label="Other", + k=best_k_w,# re-run with best k for weighted voting + #conf_threshold =best_params$threshold, use_weights = best_params$use_weights, - ignore_top = ignore_top, - verbose = verbose) + ignore_self = ignore_self, + verbose = verbose, + track_neighbors = TRUE) + + pred_with_truth_full = bind_cols(outs$df %>% select(sample_id, !!sym(truth_column)), pred) + + best_pred_w = process_votes(df=pred_with_truth_full, + group_labels=truth_classes, + k=best_k_w) + + best_pred_w = best_pred_w %>% + mutate(DLBCLone_w = ifelse(score_ratio >= best_w_purity | + top_group_score > best_w_score_thresh, + by_score, + "Other")) + + if(optimize_for_other){ + + pred = mutate(pred,predicted_label_optimized = ifelse(other_score > best_params$threshold_outgroup, + "Other", + predicted_label)) + }else{ + pred = mutate(pred,predicted_label_optimized = predicted_label) + } + + #new naming convention + pred = mutate(pred, DLBCLone_i = predicted_label, DLBCLone_io = predicted_label_optimized) + best_pred_w = rename(best_pred_w, DLBCLone_wo = DLBCLone_w) #optimized + best_pred_w = rename(best_pred_w, DLBCLone_w = by_score) #greedy + #check accuracy again + print(dim(pred)) xx_d = bind_cols(outs$df,pred) + print(dim(xx_d)) + xx_d = left_join(xx_d, select(best_pred_w, sample_id, score_ratio, top_group_score, DLBCLone_w, DLBCLone_wo), by="sample_id") + acc_check_w = report_accuracy(xx_d,pred = "DLBCLone_w") + acc_check_wo = report_accuracy(xx_d,pred = "DLBCLone_wo") + message(paste("DLBCLone_w accuracy:", + round(acc_check_w$mean_balanced_accuracy,3), + "DLBCLone_wo accuracy:", + round(acc_check_wo$mean_balanced_accuracy,3), + "DLBCLone_w Concordance:", round(acc_check_w$no_other,3), + "DLBCLone_wo Concordance:", round(acc_check_wo$no_other,3))) to_ret = list(params=results, best_params = best_params, - model=best_fit$model, - features=best_fit$features, - df=outs$df, - predictions=xx_d) + model=umap_out$model, + features=umap_out$features, + k_DLBCLone_i = best_params$k, + threshold_DLBCLone_i = best_params$threshold, + theshold_outgroup_DLBCLone_i = best_params$threshold_outgroup, + k_DLBCLone_w = best_k_w, + purity_DLBCLone_w = best_w_purity, + score_thresh_DLBCLone_w = best_w_score_thresh, + df=outs$df, + predictions=xx_d, + best_pred_w = best_pred_w, + truth_classes = truth_classes, + optimize_for_other = optimize_for_other, + truth_column = truth_column, + type = "DLBCLone_optimize_params") if(!"Other" %in% truth_classes && n_other > 0){ to_ret[["predictions_other"]] = xx_o to_ret[["predictions_combined"]] = bind_rows(xx_o,best_pred) } + #TODO: roll components of umap_out into to_ret + if(any(!names(outs) %in% names(to_ret))){ + missing_names = names(outs)[!names(outs) %in% names(to_ret)] + for(mn in missing_names){ + to_ret[[mn]] = outs[[mn]] + } + } return(to_ret) } - - #' Weighted k-nearest neighbor with confidence estimate #' #' @param train_coords Data frame containing coordinates for samples with known @@ -545,25 +2506,44 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, #' @param verbose Whether to print verbose outputs to console #' @param use_weights Set to FALSE for all neigbors to have equal weight when #' calculating the confidence -#' @param ignore_top Set to TRUE to avoid considering a nearest neighbor with -#' distance = 0. This is usually only relevant when re-classifying labeled -#' samples to estimate overall accuracy #' @param track_neighbors Set to TRUE to include details for the nearest neighbors of each sample +#' @param separate_other Set to TRUE to treat the "Other" class separately +#' when calculating the confidence. +#' @param max_neighbors Maximum number of neighbors to consider for each sample. Default 500. #' -#' @returns data frame with labels and confidence values for rows in test_coords +#' @return Data frame with rows = test samples and columns: +#' \item{predicted_label}{the predicted class} +#' \item{confidence}{predicted class weight / total weight} +#' If \code{track_neighbors = TRUE}, additional columns: +#' \item{other_score}{relative weight of outgroup vs predicted class} +#' \item{neighbor_id}{comma-separated neighbor sample IDs} +#' \item{neighbor}{comma-separated neighbor indices (in train order)} +#' \item{distance}{comma-separated neighbor distances} +#' \item{label}{comma-separated neighbor labels (in-group only if separate_other=TRUE)} +#' \item{vote_labels}{comma-separated unique labels contributing to weights} +#' \item{weighted_votes}{comma-separated weights per \code{vote_labels}} +#' \item{neighbors_other}{count of outgroup neighbors closer than the farthest in-group neighbor} +#' \item{other_weighted_votes}{sum of outgroup weights closer than the farthest in-group neighbor} +#' \item{total_w}{sum of weights for in-group neighbors used} +#' \item{pred_w}{weight supporting the predicted class} +#' +#' @import FNN dplyr tibble tidyr #' @export #' weighted_knn_predict_with_conf <- function(train_coords, train_labels, test_coords, k, - epsilon = 1e-5, + epsilon = 0.1, conf_threshold = NULL, - na_label = "Other", + other_class = "Other", verbose = FALSE, use_weights = TRUE, - ignore_top = FALSE, - track_neighbors = FALSE) { + ignore_self = TRUE, + track_neighbors = TRUE, + separate_other = TRUE, + max_neighbors = 500) { + if (nrow(train_coords)==0 || nrow(test_coords) == 0) { print("train_coords:") print(nrow(train_coords)) @@ -571,25 +2551,36 @@ weighted_knn_predict_with_conf <- function(train_coords, print(nrow(test_coords)) stop("train_coords and test_coords must be data frames with at least one row") } - - nn <- get.knnx(train_coords, test_coords, k) + # get the 100 nearest neighbors + nn <- get.knnx(train_coords, test_coords, max_neighbors) all_neighbors = data.frame() preds <- character(nrow(test_coords)) confs <- numeric(nrow(test_coords)) train_labels = as.character(train_labels) for (i in 1:nrow(test_coords)) { + + self = rownames(test_coords)[i] + + if(verbose){ + print(paste("index:",i)) + print(test_coords[i,]) + } neighbors <- nn$nn.index[i, ] distances <- nn$nn.dist[i, ] - if(ignore_top){ - #ignore a neighbor if it has identical V1 and V2 - distances <- nn$nn.dist[i, ] - if(distances[1] == 0){ - neighbors = neighbors[-1] - distances = distances[-1] - } - + if(ignore_self){ + nns = length(neighbors) + neighbor_ids = rownames(train_coords)[neighbors] + + neighbors <- neighbors[neighbor_ids != self] + distances <- distances[neighbor_ids != self] + neighbor_ids = neighbor_ids[neighbor_ids != self] + + neighbor_labels <- train_labels[neighbors] + } + + distances = distances + epsilon weights <- 1 / distances if(use_weights){ @@ -612,55 +2603,868 @@ weighted_knn_predict_with_conf <- function(train_coords, print(neighbor_labels) } # Remove NAs (just in case) + other_mask <- (neighbor_labels == other_class) + other_dists <- distances[other_mask] + other_ids <- neighbor_ids[other_mask] + other_labels <- neighbor_labels[other_mask] valid <- !is.na(neighbor_labels) + if(!all(other_labels == other_class)){ + stop("logic error: not all other_labels are other_class") + } + #num_other_neighbors = sum(neighbor_labels == "Other") + #other_dists = distances[neighbor_labels == "Other"] + if (separate_other) valid <- valid & !other_mask + neighbor_labels <- neighbor_labels[valid] weights <- weights[valid] + distances <- distances[valid] + neighbors <- neighbors[valid] + #number of neighbours should be, at least, k - 1. If less than that, warn the user + if(length(neighbors) < k-1){ + print(paste("Warning: number of neighbors is less than k-1.")) + print(paste("i:", i,"k:",k)) + print(paste("num_neighbors:",length(neighbors))) + print(table(valid)) + } + #now take the first k neighbors + if(length(neighbor_labels) > k){ + neighbor_labels = neighbor_labels[1:k] + weights = weights[1:k] + distances = distances[1:k] + neighbors = neighbors[1:k] + } + + others_closer = which(other_dists < max(distances)) + + others_distances = other_dists[others_closer] + other_ids = other_ids[others_closer] + if(use_weights){ + others_weights = 1 / others_distances + }else{ + others_weights = rep(1,length(others_closer)) + } + neighbors_other = length(others_closer) + other_weighted_votes = sum(others_weights) + mean_other_dist = mean(others_distances) + if (length(neighbor_labels) == 0) { - preds[i] <- na_label - confs[i] <- NA - next + preds[i] <- "Other" + confs[i] <- 1 + if(track_neighbors){ + + rel_other = 10 + neighbor_info <- data.frame( + other_score = rel_other, + neighbor_id = paste(rownames(train_coords)[neighbors],collapse=","), + neighbor = paste(neighbors,collapse=","), + other_neighbor = paste(other_ids,collapse=","), + distance = paste(round(distances, 3),collapse=","), + label = paste(neighbor_labels,collapse=","), + weighted_votes = "", + neighbors_other = neighbors_other, + other_weighted_votes = 0, + mean_other_dist = mean_other_dist, + total_w = 1, + pred_w = 2 + + ) + all_neighbors = bind_rows(all_neighbors, neighbor_info) + next; + } } weighted_votes <- tapply(weights, neighbor_labels, sum) - if(track_neighbors){ - # Create a data frame to store neighbors, distances, and weights - neighbor_info <- data.frame( - neighbor = paste(neighbors,collapse=","), - distance = paste(round(distances, 3),collapse=","), - #weight = paste(round(weights, 3),collapse=","), - label = paste(neighbor_labels,collapse=","), - weighted_votes = paste(weighted_votes,collapse=",") - ) - all_neighbors = bind_rows(all_neighbors, neighbor_info) - } + if (length(weighted_votes) == 0) { - preds[i] <- na_label + preds[i] <- other_class confs[i] <- NA + total_weight <- 0 + pred_weight <- 0 } else { + #print(weighted_votes) predicted_label <- names(which.max(weighted_votes)) total_weight <- sum(weighted_votes) pred_weight <- weighted_votes[predicted_label] confidence <- pred_weight / total_weight # Confidence thresholding - if (!is.null(conf_threshold) && confidence < conf_threshold) { - preds[i] <- na_label - confs[i] <- confidence - } else { - preds[i] <- predicted_label - confs[i] <- confidence + #if (!is.null(conf_threshold) && confidence < conf_threshold) { + # preds[i] <- na_label + # confs[i] <- confidence + #} else { + preds[i] <- predicted_label + confs[i] <- confidence + #} + } + + if(track_neighbors){ + # Create a data frame to store neighbors, distances, and weights + if(separate_other){ + rel_other = other_weighted_votes / pred_weight + }else{ + rel_other = 0 } + + neighbor_info <- data.frame( + other_score = rel_other, + neighbor_id = paste(rownames(train_coords)[neighbors],collapse=","), + neighbor = paste(neighbors,collapse=","), + distance = paste(round(distances, 3),collapse=","), + label = paste(neighbor_labels,collapse=","), + other_neighbor = paste(other_ids,collapse=","), + vote_labels = paste(names(weighted_votes),collapse=","), + weighted_votes = paste(weighted_votes,collapse=","), + neighbors_other = neighbors_other, + other_weighted_votes = other_weighted_votes, + total_w = total_weight, + pred_w = pred_weight + + ) + all_neighbors = bind_rows(all_neighbors, neighbor_info) } + } to_return = data.frame( - predicted_label = factor(preds), + predicted_label = preds, confidence = confs) if(track_neighbors){ - print(dim(to_return)) - print(dim(all_neighbors)) + + #check for any missing points + if(!nrow(to_return)==nrow(all_neighbors)){ + print("mismatch in row number for to_return and all_neighbors") + print(nrow(to_return)) + print(nrow(all_neighbors)) + print(head(to_return)) + print(head(all_neighbors)) + stop("") + } to_return = bind_cols(to_return,all_neighbors) + if(nrow(to_return ) != nrow(test_coords)){ + print("mismatch in row number for to_return and test_coords") + print(nrow(to_return)) + print(nrow(test_coords)) + print(head(to_return)) + print(head(test_coords)) + stop("") + } + rownames(to_return) <- rownames(test_coords) + } + return(to_return) +} + +#' @export +prepare_single_sample_DLBCLone <- function(optimized_model,seed=12345){ + if(!optimized_model$type == "DLBCLone_optimize_params"){ + stop("Input must be the output of predict_single_sample_DLBCLone") + } + if(!"projection" %in% names(optimized_model)){ + projection <- make_and_annotate_umap( + df = optimized_model$features, + umap_out = optimized_model, + ret_model = FALSE, + seed = seed, + join_column = "sample_id", + na_vals = optimized_model$best_params$na_option + ) + }else{ + stop("Model already contains a projection. Nothing to do!") + } + optimized_model$projection = projection + return(optimized_model) +} + +#' Predict class for a single sample without using umap_transform and plot result of classification +#' +#' @param test_df Data frame containing the mutation status of the test sample +#' @param train_metadata Metadata for training samples with truth labels in lymphgen column +#' @param optimized_model list of parameters from DLBCLone_optimize_params, neccessary UMAP output from a +#' previous, Data frame with the best parameters. useful for reproducibility. +#' @param truth_classes Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2") +#' @param seed Random seed for reproducibility +#' @param max_neighbors Maximum number of neighbors to consider for each sample. Default 500. +#' +#' @returns a list of data frames with the predictions, the UMAP input, the model, and a ggplot object +#' +#' @import dplyr tibble ggplot2 readr +#' @export +#' +#' @examples +#' +#' library(GAMBLR.predict) +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) +#' +#' optimize_params <- DLBCLone_optimize_params( +#' all_full_status, +#' dlbcl_meta, +#' umap_out = make_umap, +#' truth_classes = c("MCD","EZB","BN2","N1","ST2","Other"), +#' optimize_for_other = T, +#' min_k=5, +#' max_k=23 +#' ) +#' +#' predict_single <- predict_single_sample_DLBCLone( +#' test_df = optimize_params$features[1,], +#' train_metadata = dlbcl_meta, +#' optimized_model = optimize_params +#' ) +#' +predict_single_sample_DLBCLone <- function( + test_df, + train_metadata, + optimized_model = NULL, + truth_classes = c("EZB","MCD","ST2","N1","BN2","Other"), + seed = 12345, + max_neighbors = 500 +){ + if(is.null(optimized_model)){ + stop("optimized_model the output of DLBCLone_optimize_params is a required argument. Please update your code accordingly") + } + + if(any(!colnames(test_df) %in% colnames(optimized_model$features))){ + stop( + "test_df should not contain any features that are not in training data. ", + "Please check the column names of test_df and optimized_model$features." + ) + } + + bad_samps <- rowSums(test_df) == 0 + bad_test_df <- test_df[bad_samps, ] %>% + rownames_to_column("sample_id") %>% + select(sample_id) + test_df <- test_df[!bad_samps, ] + + test_id <- test_df %>% + rownames_to_column("sample_id") %>% + pull(sample_id) + + train_coords = dplyr::filter( + optimized_model$df + ) %>% + select(sample_id,V1,V2) %>% + column_to_rownames("sample_id") + + train_df_proj = dplyr::filter( + optimized_model$df + ) %>% + select(sample_id,V1,V2) %>% + left_join( #Join to the incoming metadata rather than trusting the metadata in the projection + train_metadata %>% + select( + sample_id, + lymphgen + ), + by = "sample_id" + ) + + train_labels = train_df_proj %>% + pull(lymphgen) + + #Obtain UMAP coordinates for the test sample(s) + print(paste("num rows:",nrow(test_df))) + test_projection <- make_and_annotate_umap( + df = test_df, + umap_out = optimized_model, + ret_model = FALSE, + seed = seed, + join_column = "sample_id", + na_vals = optimized_model$best_params$na_option + ) + + test_coords = dplyr::filter( + test_projection$df, + sample_id %in% test_id + ) %>% + select(sample_id,V1,V2) %>% + column_to_rownames("sample_id") + + test_pred = weighted_knn_predict_with_conf( + train_coords = train_coords, + train_labels = train_labels, + test_coords = test_coords, + k = optimized_model$best_params$k, + conf_threshold = optimized_model$best_params$threshold, + na_label = "Other", + use_weights = optimized_model$best_params$use_weights, + max_neighbors = max_neighbors + ) + + test_pred = rownames_to_column(test_pred, var = "sample_id") + if(!is.null(optimized_model)){ + test_pred = mutate( + test_pred, + DLBCLone_i = predicted_label, + DLBCLone_io = ifelse( + other_score > optimized_model$best_params$threshold_outgroup, + "Other", + predicted_label + ) + ) + + }else{ + test_pred = mutate( + test_pred, + DLBCLone_i = predicted_label, + DLBCLone_io = ifelse( + other_score > optimized_model$best_params$threshold_outgroup, + "Other", + predicted_label + ) + ) } + + anno_umap = select(test_projection$df, sample_id, V1, V2) + + anno_out = left_join(test_pred,anno_umap,by="sample_id") %>% + mutate(label = paste(sample_id,predicted_label,round(confidence,3))) + + anno_out = anno_out %>% + mutate( + V1 = as.numeric(V1), + V2 = as.numeric(V2), + label = as.character(label) + ) + + predictions_train_df = filter( + optimized_model$df + ) %>% + select(sample_id, V1, V2) + + predictions_test_df = left_join(test_pred, test_projection$df, by = "sample_id") + + predictions_df = bind_rows( + predictions_train_df %>% select(sample_id, V1, V2), + predictions_test_df %>% select(sample_id, V1, V2) + ) + + if(!is.null(optimized_model)){ + best_w_purity = optimized_model$best_w_purity + best_w_score_thresh = optimized_model$best_w_score_thresh + + predictions_test_df = process_votes( + df=predictions_test_df, + group_labels=optimized_model$truth_classes, + k=optimized_model$k_DLBCLone_w + ) + + predictions_test_df = predictions_test_df %>% + mutate( + DLBCLone_w = by_score, + DLBCLone_wo = ifelse( + score_ratio >= optimized_model$purity_DLBCLone_w | top_group_score > optimized_model$score_thresh_DLBCLone_w, + by_score, + "Other" + ) + ) + } + + predictions_test_df = predictions_test_df %>% + bind_rows( + bad_test_df %>% + mutate( + predicted_label = "Other", + confidence = 1 + ) + ) + + combined_df <- bind_rows( + optimized_model$features %>% rownames_to_column("sample_id"), + test_df %>% rownames_to_column("sample_id") + ) %>% + distinct(sample_id, .keep_all = TRUE) %>% # keep first occurrence + column_to_rownames("sample_id") + + if(any(!colnames(optimized_model$features) %in% colnames(test_df))){ + message("filling in missing features in test_df with zeros") + combined_df[is.na(combined_df)] = 0 + } + + to_return = list( + predictions = predictions_test_df, + df = predictions_df, + anno_df = predictions_df %>% left_join(.,train_metadata,by="sample_id"), + anno_out = anno_out, + features_df = combined_df, + type = "predict_single_sample_DLBCLone", + optimized_predictions = optimized_model$predictions + ) + return(to_return) } +#' model storage for DLBCLone outputs +#' +#' @param optimized_params List containing the optimized parameters from DLBCLone_optimize_params +#' @param path Path to save the files +#' @param name_prefix Prefix for the saved files, all files will be in path and start with name_prefix +#' +#' @returns saves the files to the specified path +#' +#' @import uwot readr +#' +#' @export +#' +#' @examples +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' DLBCLone_save_optimized( +#' optimzied_params = optimized_params, +#' path="/save_optimized/trial_folder", +#' name_prefix="test_A" +#' ) +#' } +#' +DLBCLone_save_optimized = function( + optimized_params=NULL, + path="models/", + name_prefix="test" +){ + # all files will be in path and start with name_prefix + prefix = paste0(path,"/",name_prefix) + + out_mut = paste0(prefix,"_mutation_status_df.rds") + saveRDS(optimized_params$features,file=out_mut) + + out_meta = paste0(prefix,"_metadata.tsv") + write_tsv(optimized_params$df,file=out_meta) + + out_param = paste0(prefix,"_optimized_best_params.rds") + saveRDS(optimized_params$best_params,file=out_param) + + out_model = paste0(prefix,"_optimized_uwot.rds") + save_uwot(optimized_params$model,file=out_model) + + out_pred = paste0(prefix,"_optimized_pred.tsv") + write_tsv(optimized_params$predictions,file=out_pred) + + out_classes = paste0(prefix,"_classes.rds") + saveRDS(optimized_params$truth_classes,file=out_classes) + + out_k_w = paste0(prefix,"_optimized_k_w.rds") + saveRDS(optimized_params$k_DLBCLone_w,file=out_k_w) + + purity_w = paste0(prefix,"_purity_w.rds") + saveRDS(optimized_params$purity_DLBCLone_w,file=purity_w) + + score_thresh_w = paste0(prefix,"_score_thresh_w.rds") + saveRDS(optimized_params$score_thresh_DLBCLone_w,file=score_thresh_w) +} + +#' load previously saved DLBCLone model and parameters +#' +#' @param path Path to open saved the files +#' @param name_prefix Prefix of the saved files, all files will be in path and start with name_prefix +#' +#' @returns saves the files to the specified path +#' +#' @import uwot readr dplyr +#' @export +#' +#' @examples +#' +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' load_optimized <- DLBCLone_load_optimized( +#' path="/save_optimized/trial_folder", +#' name_prefix="test_A" +#' ) +#' } +#' +DLBCLone_load_optimized <- function( + path="models/", + name_prefix="test" +){ + #all files will be in path and start with name_prefix + prefix = paste0(path,"/",name_prefix) + + load_mut = paste0(prefix,"_mutation_status_df.rds") + load_meta = paste0(prefix,"_metadata.tsv") + load_classes = paste0(prefix,"_classes.rds") + load_param = paste0(prefix,"_optimized_best_params.rds") + load_model = paste0(prefix,"_optimized_uwot.rds") + load_pred = paste0(prefix,"_optimized_pred.tsv") + load_k_w = paste0(prefix,"_optimized_k_w.rds") + load_purity_w = paste0(prefix,"_purity_w.rds") + load_score_thresh_w = paste0(prefix,"_score_thresh_w.rds") + + required_files <- c(load_mut, load_meta, load_classes, load_param, load_model, load_pred, load_k_w, load_purity_w, load_score_thresh_w) + + # Check existence + missing_files <- required_files[!file.exists(required_files)] + if (length(missing_files) > 0) { + stop("The following required files are missing:\n", paste(missing_files, collapse = "\n")) + } + + mut_df <- readRDS(load_mut) + metadata <- read_tsv(load_meta) %>% + mutate(lymphgen = as.factor(lymphgen)) + classes <- readRDS(load_classes) + best_params <- readRDS(load_param) + uwot_model <- load_uwot(load_model) + predictions <- read_tsv(load_pred) + k_DLBCLone_w <- readRDS(load_k_w) + best_w_purity <- readRDS(load_purity_w) + best_w_score_thresh <- readRDS(load_score_thresh_w) + + return(list( + features = mut_df, + df = metadata, + truth_classes = classes, + best_params = best_params, + model = uwot_model, + predictions = predictions, + k_DLBCLone_w = k_DLBCLone_w, + purity_DLBCLone_w = best_w_purity, + score_thresh_DLBCLone_w = best_w_score_thresh + )) +} + +#' Train a Gaussian Mixture Model for DLBCLone Classification +#' +#' Fits a supervised Gaussian mixture model (GMM) to UMAP-projected data for DLBCLone subtypes, excluding samples labeled "Other". +#' Assigns class predictions and optionally reclassifies samples as "Other" based on probability and density thresholds. +#' +#' @param umap_out List. Output from \code{make_and_annotate_umap}, containing a data frame with UMAP coordinates and truth labels. +#' @param probability_threshold Numeric. Minimum posterior probability required to assign a class (default: 0.5). +#' @param density_max_threshold Numeric. Minimum maximum density required to assign a class (default: 0.05). +#' @param truth_column Name for the original class column in the predictions data frame (default: "lymphgen"). +#' @param cohort Optional character. Cohort label to annotate predictions. +#' @param truth_classes Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2","Other"). +#' +#' @details +#' - Uses \code{MclustDA} to fit a supervised mixture model to the UMAP coordinates (V1, V2) and class labels. +#' - Predicts class membership and computes per-class densities for each sample. +#' - Samples with low maximum probability or density are reclassified as "Other". +#' - Returns both raw and thresholded class assignments, respectively under the columns DLBCLone_g and DLBCLone_go. +#' +#' @return A list with: +#' \item{gaussian_mixture_model}{Fitted \code{MclustDA} model object} +#' \item{predictions}{Data frame with sample IDs, UMAP coordinates, true labels, predicted classes, and thresholded assignments} +#' \item{probability_threshold}{Probability threshold used for "Other" assignment} +#' +#' @import mclust dplyr tibble +#' @export +#' +#' @examples +#' +#' library(GAMBLR.predict) +#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) +#' +#' result <- DLBCLone_train_mixture_model(umap_out = make_umap) +#' head(result$predictions) +#' +DLBCLone_train_mixture_model = function(umap_out, + probability_threshold = 0.5, + density_max_threshold = 0.05, + truth_column = "lymphgen", + cohort = NULL, + truth_classes = c("EZB","MCD","ST2","N1","BN2","Other") + ){ + + df = umap_out$df %>% select(sample_id,!!sym(truth_column),V1,V2) %>% filter(!!sym(truth_column) != "Other") + df = filter(df,!!sym(truth_column) %in% truth_classes) + #mont_test_proj = montreal_gambl_c_mu$df %>% select(sample_id,lymphgen,V1,V2) + + df[[truth_column]] <- as.factor(df[[truth_column]]) + + #model without Others + gmm_supervised <- MclustDA( + df[, c("V1", "V2")], + class = df[[truth_column]], + modelType = "MclustDA" # instead of "EDDA" + ) +df = umap_out$df %>% select(sample_id,!!sym(truth_column),V1,V2) #%>% filter(!!sym(truth_column) != "Other") + +pred <- predict(gmm_supervised, newdata = df[, c("V1", "V2")]) + +densities_list <- lapply(names(gmm_supervised$models), function(class_name) { + model <- gmm_supervised$models[[class_name]] + + modelName <- model$modelName + params <- model$parameters + + # Compute log-densities for each component + logdens <- cdens( + data = df[, c("V1", "V2")], + modelName = modelName, + parameters = params, + logarithm = TRUE + ) + + # cdens can return a matrix (samples x components); sum components + class_density <- exp(logdens) # back to raw densities + if (is.matrix(class_density)) { + class_density <- rowSums(class_density) + } + + class_density +}) + +densities <- do.call(cbind, densities_list) +colnames(densities) <- names(gmm_supervised$models) + +max_density <- apply(densities, 1, max) +max_prob <- apply(pred$z, 1, max) + +df$DLBCLone_g = pred$classification + +df$DLBCLone_go <- ifelse( + (max_prob < probability_threshold) | (max_density < density_max_threshold), + "Other", + as.character(pred$classification) +) + +df$cohort = cohort +to_return = list(gaussian_mixture_model = gmm_supervised, + predictions = df, + truth_classes = truth_classes, + truth_column = truth_column, + probability_threshold = probability_threshold + ) + +return(to_return) +} + +#' Predict DLBCLone Class Membership Using a Trained Gaussian Mixture Model +#' +#' Applies a previously trained supervised Gaussian mixture model (GMM) to UMAP-projected data for DLBCLone subtypes. +#' Assigns class predictions and optionally reclassifies samples as "Other" based on probability and density thresholds. +#' +#' @param model Fitted \code{MclustDA} model object, as returned by \code{DLBCLone_train_mixture_model}. +#' @param umap_out List. Output from \code{make_and_annotate_umap}, containing a data frame with UMAP +#' coordinates for the samples to be classified with the model. This must be projected using the same UMAP model +#' that was generated using the training data. +#' @param probability_threshold Numeric. Minimum posterior probability required to assign a class (default: 0.5). +#' @param density_max_threshold Numeric. Minimum maximum density required to assign a class (default: 0.05). +#' @param cohort Optional character. Cohort label to annotate predictions. +#' +#' @details +#' - Uses the provided \code{MclustDA} model to predict class membership for each sample in the UMAP projection. +#' - Computes per-class densities and posterior probabilities for each sample. +#' - Samples with low maximum probability or density are reclassified as "Other". +#' - Returns both raw and thresholded class assignments, respectively under the columns DLBCLone_g and DLBCLone_go. +#' +#' @return A list with: +#' \item{gaussian_mixture_model}{Fitted \code{MclustDA} model object} +#' \item{predictions}{Data frame with sample IDs, UMAP coordinates, predicted classes, and thresholded assignments} +#' \item{probability_threshold}{Probability threshold used for "Other" assignment} +#' +#' @import mclust dplyr tibble +#' @export +#' +#' @examples +#' # Predict on new UMAP data using a trained mixture model: +#' +#' library(GAMBLR.predict) +#' +#' all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) %>% +#' tibble::column_to_rownames("sample_id") +#' +#' dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) +#' +#' mixture_model <- DLBCLone_train_mixture_model(umap_out = make_umap) +#' +#' result <- DLBCLone_predict_mixture_model( +#' model=mixture_model, +#' umap_out=make_umap +#' ) +#' head(result$predictions) +#' +#' +DLBCLone_predict_mixture_model = function(model, + umap_out, + probability_threshold = 0.5, + density_max_threshold = 0.05, + cohort = NULL + ){ + df = umap_out$df %>% select(sample_id,V1,V2) + gmm_supervised = model +pred <- predict(gmm_supervised, newdata = df[, c("V1", "V2")]) + +densities_list <- lapply(names(gmm_supervised$models), function(class_name) { + model <- gmm_supervised$models[[class_name]] + modelName <- model$modelName + params <- model$parameters + + # Compute log-densities for each component + logdens <- cdens( + data = df[, c("V1", "V2")], + modelName = modelName, + parameters = params, + logarithm = TRUE + ) + + # cdens can return a matrix (samples x components); sum components + class_density <- exp(logdens) # back to raw densities + if (is.matrix(class_density)) { + class_density <- rowSums(class_density) + } + + class_density +}) + +densities <- do.call(cbind, densities_list) +colnames(densities) <- names(gmm_supervised$models) + +max_density <- apply(densities, 1, max) +max_prob <- apply(pred$z, 1, max) + + +df$DLBCLone_g = pred$classification + + +df$DLBCLone_go <- ifelse( + (max_prob < probability_threshold) | (max_density < density_max_threshold), + "Other", + as.character(pred$classification) +) + + +df$cohort = cohort +to_return = list(gaussian_mixture_model = gmm_supervised, + predictions = df, + probability_threshold = probability_threshold + ) +return(to_return) +} + +# Some comment here +# + +old_process_votes <- function(df, + raw_col = "label", + group_labels = c("EZB", "MCD", "ST2", "BN2", "N1", "Other"), + vote_labels_col = "vote_labels", + k, + other_vote_multiplier = 2, + score_purity_requirement = 1, + weighted_votes_col = "weighted_votes") { + if(missing(k)){ + stop("k value is required") + } + score_thresh = 2 * k + + count_labels_in_string <- function(string, labels) { + tokens <- str_split(string, ",")[[1]] + map_int(labels, ~ sum(tokens == .x)) + } + + extract_weighted_scores <- function(label_str, vote_str, labels) { + lbls <- str_split(label_str, ",")[[1]] + votes <- as.numeric(str_split(vote_str, ",")[[1]]) + map_dbl(labels, ~ sum(votes[lbls == .x])) %>% + set_names(paste0(labels, "_score")) + } + + get_top_score_group <- function(label_str, vote_str, labels) { + lbls <- str_split(label_str, ",")[[1]] + votes <- as.numeric(str_split(vote_str, ",")[[1]]) + scores_by_label <- set_names(map_dbl(labels, ~ sum(votes[lbls == .x])), labels) + top <- names(scores_by_label)[which.max(scores_by_label)] + value <- scores_by_label[[top]] + list(top_score_group = top, top_group_score = value) + } + + df_out <- df %>% + mutate(.id = row_number()) %>% + rowwise() %>% + mutate( + # 1) counts + counts = list( + set_names( + count_labels_in_string(.data[[raw_col]], group_labels), + paste0(group_labels, "_NN_count") + ) + ), + top_group = { + cnts <- count_labels_in_string(.data[[raw_col]], group_labels) + group_labels[which.max(cnts)] + }, + + # 2) scores + scores = list( + if (!is.null(vote_labels_col) && !is.null(weighted_votes_col)) { + extract_weighted_scores( + .data[[vote_labels_col]], + .data[[weighted_votes_col]], + group_labels + ) + } else { + set_names(rep(0, length(group_labels)), paste0(group_labels, "_score")) + } + ), + + # 3) top score group summary + score_summary = list( + if (!is.null(vote_labels_col) && !is.null(weighted_votes_col)) { + get_top_score_group( + .data[[vote_labels_col]], + .data[[weighted_votes_col]], + group_labels + ) + } else { + list(top_score_group = NA_character_, top_group_score = NA_real_) + } + ) + ) %>% + ungroup() %>% + unnest_wider(counts) %>% + unnest_wider(scores) %>% + unnest_wider(score_summary) %>% + rowwise() %>% +mutate( + top_group_count = get(paste0(top_group, "_NN_count")) +) %>% +ungroup() + + # Optional adjustments based on external columns (if present) + if ("neighbors_other" %in% colnames(df)) { + df_out <- df_out %>% + mutate(Other_count = neighbors_other) + } + + if ("other_weighted_votes" %in% colnames(df)) { + df_out <- df_out %>% + mutate(Other_score = other_weighted_votes) + } + # Optional override of top group if Other dominates by a fold + if (all(c("top_group_count", "Other_count") %in% colnames(df_out))) { + df_out <- df_out %>% + mutate(by_vote = top_group_count) %>% + mutate(by_vote_opt = ifelse(top_group_count * other_vote_multiplier > Other_count, top_group, "Other")) + } + + df_out = mutate(df_out, + by_score = top_score_group, + score_ratio = top_group_score / Other_score, + by_score_opt=ifelse(score_ratio > score_purity_requirement | top_group_score > score_thresh,top_score_group,"Other")) + + return(df_out) +} diff --git a/man/DLBCLone_KNN.Rd b/man/DLBCLone_KNN.Rd new file mode 100644 index 0000000..76c2048 --- /dev/null +++ b/man/DLBCLone_KNN.Rd @@ -0,0 +1,106 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{DLBCLone_KNN} +\alias{DLBCLone_KNN} +\title{Run DLBCLone KNN Classification} +\usage{ +DLBCLone_KNN( + features_df, + metadata, + core_features = NULL, + core_feature_multiplier = 1.5, + hidden_features = NULL, + min_k = 5, + max_k = 60, + truth_column = "lymphgen", + truth_classes = c("EZB", "BN2", "ST2", "MCD", "N1", "Other"), + other_class = "Other", + optimize_for_other = TRUE, + predict_unlabeled = FALSE, + plot_samples = NULL, + DLBCLone_KNN_out = NULL, + seed = 12345, + epsilon = 0.001, + weighted_votes = TRUE, + skip_umap = FALSE +) +} +\arguments{ +\item{features_df}{Numeric matrix/data.frame (rows = samples, cols = features). +Row names must be sample IDs.} + +\item{metadata}{Data frame with at least \code{sample_id} and the ground-truth +label column given in \code{truth_column}.} + +\item{core_features}{Character vector of feature names to upweight (optional).} + +\item{core_feature_multiplier}{Numeric multiplier for \code{core_features}.} + +\item{hidden_features}{Character vector of feature names to drop (optional).} + +\item{min_k, max_k}{Integer K range to explore when optimizing.} + +\item{truth_column}{Name of metadata column with ground-truth class labels.} + +\item{truth_classes}{Character vector of all classes to consider (including +\code{other_class} if you intend to optimize for it).} + +\item{other_class}{Name of the explicit outgroup class (default: "Other").} + +\item{optimize_for_other}{Logical; if TRUE, computes a separate "other" +score (ratio) and searches a purity threshold; if FALSE, treats all +classes symmetrically.} + +\item{predict_unlabeled}{If TRUE, re-runs KNN to classify samples that were +present in \code{features_df} but not in \code{metadata}.} + +\item{plot_samples}{Optional vector of sample_ids to keep in example plots.} + +\item{DLBCLone_KNN_out}{Optional prior result; if supplied, its learned +parameters are reused (skip optimization).} + +\item{seed}{Random seed.} + +\item{epsilon}{Small value added to distances before weighting.} + +\item{weighted_votes}{If FALSE, neighbors are unweighted (equal votes).} + +\item{skip_umap}{If TRUE, skip layout optimization plots at the end.} +} +\value{ +A list with fields including: +\item{predictions}{Per-sample vote/score summary and predicted labels} +\item{DLBCLone_k_best_k}{Best K found} +\item{DLBCLone_k_purity_threshold}{Best purity threshold (if applicable)} +\item{DLBCLone_k_accuracy}{Best accuracy metric achieved} +\item{truth_classes, truth_column}{Echoed arguments} +\item{unlabeled_predictions}{Predictions for unlabeled samples (if requested)} +\item{df}{Annotated layout for plotting (when built in this run)} +\item{plot_truth, plot_predicted}{ggplots when built in this run} +} +\description{ +Weighted KNN on a feature (mutation) matrix with optional upweighting of +user-specified "core" features, optional exclusion of "hidden" features, +and optional optimization of an explicit outgroup (e.g. "Other"). +} +\details{ +This version removes hard-coded LymphGen class names and instead derives the +in-group classes and the outgroup column name from the arguments +\code{truth_classes} and \code{other_class}. It keeps backward compatibility +for the default LymphGen-like usage. +} +\examples{ + +library(GAMBLR.predict) + +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +dlbcl_knn <- DLBCLone_KNN( + features_df = all_full_status, + metadata = dlbcl_meta +) + +} diff --git a/man/DLBCLone_KNN_predict.Rd b/man/DLBCLone_KNN_predict.Rd new file mode 100644 index 0000000..1f7b9a6 --- /dev/null +++ b/man/DLBCLone_KNN_predict.Rd @@ -0,0 +1,77 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{DLBCLone_KNN_predict} +\alias{DLBCLone_KNN_predict} +\title{Predict DLBCLone Classes for New Samples Using a Trained KNN Model} +\usage{ +DLBCLone_KNN_predict( + train_df, + test_df, + metadata, + DLBCLone_KNN_out, + mode = "batch", + truth_column = "lymphgen", + other_class = "Other" +) +} +\arguments{ +\item{train_df}{Data frame or matrix of features for training samples (rows = samples, columns = features).} + +\item{test_df}{Data frame or matrix of features for test samples to be classified.} + +\item{metadata}{Data frame with metadata for all samples, including at least a \code{sample_id} column.} + +\item{DLBCLone_KNN_out}{List. Output from a previous call to \code{DLBCLone_KNN} containing optimized parameters. (Required)} + +\item{mode}{Character. If \code{"iterative"}, runs KNN prediction for each test sample individually (recommended for stability).} + +\item{core_features}{Optional character vector of feature names to upweight in the KNN calculation.} + +\item{core_feature_multiplier}{Numeric. Multiplier to apply to core features (default: 1.5).} + +\item{hidden_features}{Optional character vector of feature names to exclude from the analysis.} +} +\value{ +A list containing the KNN prediction results for the test samples, including predicted class labels and scores. +} +\description{ +Applies a previously optimized DLBCLone KNN model to predict class labels for new (test) samples. +This function combines the training and test feature matrices, ensures feature compatibility, and uses the +parameters from a DLBCLone KNN optimization run to classify the test samples. Optionally, runs in iterative mode +for more stable results when predicting multiple samples. +} +\details{ +\itemize{ +\item Ensures that the feature columns in \code{train_df} and \code{test_df} are compatible. +\item If \code{mode = "iterative"}, runs KNN prediction for each test sample one at a time. +\item Uses the parameters (e.g., k, feature weights) from the provided \code{DLBCLone_KNN_out} object. +\item Returns the same structure as \code{DLBCLone_KNN}, with predictions for the test samples. +} +} +\examples{ +# Assuming you have run DLBCLone_KNN to get optimized parameters: +# model_out <- DLBCLone_KNN(train_features, train_metadata, ...) +# Predict on new samples: + +library(GAMBLR.predict) + +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +dlbcl_knn <- DLBCLone_KNN( + features_df = all_full_status, + metadata = dlbcl_meta +) + +predictions <- DLBCLone_KNN_predict( + train_df = all_full_status, + test_df = lyseq_validation_data, + metadata = dlbcl_meta, + DLBCLone_KNN_out = dlbcl_knn, + mode = "batch" +) + + +} diff --git a/man/DLBCLone_load_optimized.Rd b/man/DLBCLone_load_optimized.Rd new file mode 100644 index 0000000..084efdf --- /dev/null +++ b/man/DLBCLone_load_optimized.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{DLBCLone_load_optimized} +\alias{DLBCLone_load_optimized} +\title{load previously saved DLBCLone model and parameters} +\usage{ +DLBCLone_load_optimized(path = "models/", name_prefix = "test") +} +\arguments{ +\item{path}{Path to open saved the files} + +\item{name_prefix}{Prefix of the saved files, all files will be in path and start with name_prefix} +} +\value{ +saves the files to the specified path +} +\description{ +load previously saved DLBCLone model and parameters +} +\examples{ + +\dontrun{ +library(GAMBLR.predict) + +load_optimized <- DLBCLone_load_optimized( + path="/save_optimized/trial_folder", + name_prefix="test_A" +) +} + +} diff --git a/man/DLBCLone_optimize_params.Rd b/man/DLBCLone_optimize_params.Rd index 95599bf..9cb9ac3 100644 --- a/man/DLBCLone_optimize_params.Rd +++ b/man/DLBCLone_optimize_params.Rd @@ -9,45 +9,114 @@ DLBCLone_optimize_params( metadata_df, umap_out, truth_classes = c("EZB", "MCD", "ST2", "N1", "BN2", "Other"), - eval_group = "Lacy", + truth_column = "lymphgen", + optimize_for_other = FALSE, + eval_group = NULL, min_k = 3, - max_k = 30, + max_k = 23, verbose = FALSE, - seed = 12345 + seed = 12345, + maximize = "balanced_accuracy", + cap_classification_rate = 0.9, + exclude_other_for_accuracy = FALSE, + weights_opt = c(TRUE) ) } \arguments{ -\item{combined_mutation_status_df}{Data frame with one row per sample and one column per mutation} +\item{combined_mutation_status_df}{Data frame with one row per sample and +one column per mutation} -\item{metadata_df}{Data frame of metadata with one row per sample and three required columns: sample_id, dataset and lymphgen} +\item{metadata_df}{Data frame of metadata with one row per sample and +three required columns: sample_id, dataset and lymphgen} -\item{umap_out}{additional argument} +\item{umap_out}{The output of a previous run of make_and_annotate_umap. +If provided, the function will use this model to project the data +instead of re-running UMAP.} -\item{truth_classes}{Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2","Other")} +\item{truth_classes}{Vector of classes to use for training and testing. +Default: c("EZB","MCD","ST2","N1","BN2","Other")} -\item{eval_group}{Specify whether certain rows will be evaluated and held out from training rather than using all samples.} +\item{truth_column}{The column in the metadata_df that contains the true class labels. +Default: "lymphgen".} -\item{min_k}{additional argument} +\item{optimize_for_other}{Set to TRUE to optimize the threshold for +classifying samples as "Other" based on the relative proportion of +samples near the sample in UMAP space with the "Other" label. Rather than +treating Other as just another class, this will optimize the threshold for +a separate score that considers how many Other and non-Other samples are +in the neighborhood of the sample in question. This parameter will NOT change +the value in predicted_label. Instead, the predicted_label_optimized column +will contain the optimized label. Default: FALSE} -\item{max_k}{additional argument} +\item{eval_group}{If desired, use this to specify which rows will be +evaluated and held out from training rather than using all samples. +NOTE: this parameter will probably become deprecated!} -\item{verbose}{additional argument} +\item{min_k}{Starting k for knn (Default: 3)} -\item{seed}{additional argument} +\item{max_k}{Ending k for knn (Default: 33)} + +\item{verbose}{Whether to print verbose outputs to console} + +\item{seed}{Random seed to use for reproducibility (default: 12345)} + +\item{maximize}{Metric to use for optimization. Either "sensitivity" +(average sensitivity across all classes), "accuracy" +(actual accuracy across all samples) or "balanced_accuracy" (the mean of the +balanced accuracy values across all classes). Default: "balanced_accuracy"} + +\item{exclude_other_for_accuracy}{Set to TRUE to exclude the +"Other" class from the 'lymphgen' column when calculating accuracy metrics +(passed to DLBCLone_optimize_params). Default: FALSE} + +\item{weights_opt}{Vector of TRUE/FALSE values to indicate whether to use +weighted knn or not. Default: c(TRUE).} } \value{ List of data frames with the results of the parameter optimization -including the best model, the associated knn parameters and the annotated UMAP output +including the best model, the associated knn parameters and the annotated +UMAP output as a data frame. The list also includes the predictions for the +"Other" class if it was included in the training and testing. } \description{ Optimize parameters for classifying samples using UMAP and k-nearest neighbor } \examples{ -lymphgen_A53_DLBCLone = DLBCLone_optimize_params( - lgen_feat_status, #our binary feature matrix - a53_meta, #our metadata - umap_out = lymphgen_A53_all_feat_gambl, # force use existing UMAP fit - eval_group = NULL, # use all samples for evaluating accuracy - truth_classes = c("MCD","EZB","BN2","ST2","N1","A53","Other")) +library(GAMBLR.predict) + +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) + +# Aim to maximize classification of samples into non-Other class +optimize_params_F <- DLBCLone_optimize_params( + all_full_status, + dlbcl_meta, + umap_out = make_umap, + truth_classes = c("MCD","EZB","BN2","N1","ST2","Other"), + optimize_for_other = F, + min_k=5, + max_k=23 +) + +# Aim to maximize balanced accuracy while allowing samples to be +# unclassified (assigned to "Other") + +optimize_params_T <- DLBCLone_optimize_params( + all_full_status, + dlbcl_meta, + umap_out = make_umap, + truth_classes = c("MCD","EZB","BN2","N1","ST2","Other"), + optimize_for_other = T, + min_k=5, + max_k=23 +) + } diff --git a/man/DLBCLone_predict_mixture_model.Rd b/man/DLBCLone_predict_mixture_model.Rd new file mode 100644 index 0000000..b3b2db3 --- /dev/null +++ b/man/DLBCLone_predict_mixture_model.Rd @@ -0,0 +1,70 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{DLBCLone_predict_mixture_model} +\alias{DLBCLone_predict_mixture_model} +\title{Predict DLBCLone Class Membership Using a Trained Gaussian Mixture Model} +\usage{ +DLBCLone_predict_mixture_model( + model, + umap_out, + probability_threshold = 0.5, + density_max_threshold = 0.05, + cohort = NULL +) +} +\arguments{ +\item{model}{Fitted \code{MclustDA} model object, as returned by \code{DLBCLone_train_mixture_model}.} + +\item{umap_out}{List. Output from \code{make_and_annotate_umap}, containing a data frame with UMAP +coordinates for the samples to be classified with the model. This must be projected using the same UMAP model +that was generated using the training data.} + +\item{probability_threshold}{Numeric. Minimum posterior probability required to assign a class (default: 0.5).} + +\item{density_max_threshold}{Numeric. Minimum maximum density required to assign a class (default: 0.05).} + +\item{cohort}{Optional character. Cohort label to annotate predictions.} +} +\value{ +A list with: +\item{gaussian_mixture_model}{Fitted \code{MclustDA} model object} +\item{predictions}{Data frame with sample IDs, UMAP coordinates, predicted classes, and thresholded assignments} +\item{probability_threshold}{Probability threshold used for "Other" assignment} +} +\description{ +Applies a previously trained supervised Gaussian mixture model (GMM) to UMAP-projected data for DLBCLone subtypes. +Assigns class predictions and optionally reclassifies samples as "Other" based on probability and density thresholds. +} +\details{ +\itemize{ +\item Uses the provided \code{MclustDA} model to predict class membership for each sample in the UMAP projection. +\item Computes per-class densities and posterior probabilities for each sample. +\item Samples with low maximum probability or density are reclassified as "Other". +\item Returns both raw and thresholded class assignments, respectively under the columns DLBCLone_g and DLBCLone_go. +} +} +\examples{ +# Predict on new UMAP data using a trained mixture model: + +library(GAMBLR.predict) + +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) + +mixture_model <- DLBCLone_train_mixture_model(umap_out = make_umap) + +result <- DLBCLone_predict_mixture_model( + model=mixture_model, + umap_out=make_umap +) +head(result$predictions) + + +} diff --git a/man/DLBCLone_save_optimized.Rd b/man/DLBCLone_save_optimized.Rd new file mode 100644 index 0000000..ef50b4d --- /dev/null +++ b/man/DLBCLone_save_optimized.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{DLBCLone_save_optimized} +\alias{DLBCLone_save_optimized} +\title{model storage for DLBCLone outputs} +\usage{ +DLBCLone_save_optimized( + optimized_params = NULL, + path = "models/", + name_prefix = "test" +) +} +\arguments{ +\item{optimized_params}{List containing the optimized parameters from DLBCLone_optimize_params} + +\item{path}{Path to save the files} + +\item{name_prefix}{Prefix for the saved files, all files will be in path and start with name_prefix} +} +\value{ +saves the files to the specified path +} +\description{ +model storage for DLBCLone outputs +} +\examples{ + +\dontrun{ +library(GAMBLR.predict) + +DLBCLone_save_optimized( + optimzied_params = optimized_params, + path="/save_optimized/trial_folder", + name_prefix="test_A" +) +} + +} diff --git a/man/DLBCLone_summarize_model.Rd b/man/DLBCLone_summarize_model.Rd new file mode 100644 index 0000000..d9566d0 --- /dev/null +++ b/man/DLBCLone_summarize_model.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{DLBCLone_summarize_model} +\alias{DLBCLone_summarize_model} +\title{Summarize and Export DLBCLone Model Results} +\usage{ +DLBCLone_summarize_model(base_name, optimized_model) +} +\arguments{ +\item{base_name}{Character. The base name (and directory) for saving output files.} + +\item{optimized_model}{List. The output from DLBCLone optimization, containing predictions, features, and metadata.} +} +\value{ +No return value. Side effect: writes multiple PDF files to disk. +} +\description{ +Generates and saves a set of summary plots and tables for a DLBCLone model, including UMAP scatterplots, alluvial plots, and oncoplots. +Results are saved as PDF files in a directory named after the provided base name. +} +\details{ +\itemize{ +\item Creates a directory for results if it does not exist. +\item Saves UMAP scatterplots for all samples and for non-"Other" samples. +\item Generates alluvial plots for different DLBCLone predictions +\item Exports an oncoplot summarizing mutation and classification results. +\item Uses \code{make_umap_scatterplot}, \code{make_alluvial}, and \code{prettyOncoplot} for visualization. +} +} +\examples{ +\dontrun{ +DLBCLone_summarize_model("Full_geneset_unweighted", optimized_model) +} + +} diff --git a/man/DLBCLone_train_mixture_model.Rd b/man/DLBCLone_train_mixture_model.Rd new file mode 100644 index 0000000..d103272 --- /dev/null +++ b/man/DLBCLone_train_mixture_model.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{DLBCLone_train_mixture_model} +\alias{DLBCLone_train_mixture_model} +\title{Train a Gaussian Mixture Model for DLBCLone Classification} +\usage{ +DLBCLone_train_mixture_model( + umap_out, + probability_threshold = 0.5, + density_max_threshold = 0.05, + truth_column = "lymphgen", + cohort = NULL, + truth_classes = c("EZB", "MCD", "ST2", "N1", "BN2", "Other") +) +} +\arguments{ +\item{umap_out}{List. Output from \code{make_and_annotate_umap}, containing a data frame with UMAP coordinates and truth labels.} + +\item{probability_threshold}{Numeric. Minimum posterior probability required to assign a class (default: 0.5).} + +\item{density_max_threshold}{Numeric. Minimum maximum density required to assign a class (default: 0.05).} + +\item{truth_column}{Name for the original class column in the predictions data frame (default: "lymphgen").} + +\item{cohort}{Optional character. Cohort label to annotate predictions.} + +\item{truth_classes}{Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2","Other").} +} +\value{ +A list with: +\item{gaussian_mixture_model}{Fitted \code{MclustDA} model object} +\item{predictions}{Data frame with sample IDs, UMAP coordinates, true labels, predicted classes, and thresholded assignments} +\item{probability_threshold}{Probability threshold used for "Other" assignment} +} +\description{ +Fits a supervised Gaussian mixture model (GMM) to UMAP-projected data for DLBCLone subtypes, excluding samples labeled "Other". +Assigns class predictions and optionally reclassifies samples as "Other" based on probability and density thresholds. +} +\details{ +\itemize{ +\item Uses \code{MclustDA} to fit a supervised mixture model to the UMAP coordinates (V1, V2) and class labels. +\item Predicts class membership and computes per-class densities for each sample. +\item Samples with low maximum probability or density are reclassified as "Other". +\item Returns both raw and thresholded class assignments, respectively under the columns DLBCLone_g and DLBCLone_go. +} +} +\examples{ + +library(GAMBLR.predict) + +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) + +result <- DLBCLone_train_mixture_model(umap_out = make_umap) +head(result$predictions) + +} diff --git a/man/DLBCLone_train_test_plot.Rd b/man/DLBCLone_train_test_plot.Rd index b394ffe..b28d581 100644 --- a/man/DLBCLone_train_test_plot.Rd +++ b/man/DLBCLone_train_test_plot.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/umap.R \name{DLBCLone_train_test_plot} \alias{DLBCLone_train_test_plot} -\title{Plot the result of a lymphgen classification} +\title{Plot the result of a DLBCLone classification} \usage{ DLBCLone_train_test_plot( test_df, @@ -13,9 +13,10 @@ DLBCLone_train_test_plot( annotate_accuracy = FALSE, classes = c("BN2", "ST2", "MCD", "EZB", "N1"), label_offset = 2, - title1 = "GAMBL", - title2 = "predicted_class_for_HighConf", - title3 = "predicted_class_for_Other" + title1 = "Original Class", + title2 = "DLBCLone Predicted Class", + title3 = "DLBCLone Predicted Class (Other)", + base_size = 1 ) } \arguments{ @@ -45,10 +46,14 @@ DLBCLone_train_test_plot( a ggplot object } \description{ -Plot the result of a lymphgen classification +Plot the result of a DLBCLone classification } \examples{ #add the dataset name to the metadata if it's not already there (required for the plot work) + +\dontrun{ +library(GAMBLR.predict) + lymphgen_A53_DLBCLone$df$dataset = "GAMBL" DLBCLone_train_test_plot( @@ -59,5 +64,6 @@ DLBCLone_train_test_plot( details = lymphgen_A53_DLBCLone$best_params, classes = c("MCD","EZB","BN2","ST2","N1","A53","Other"), annotate_accuracy=TRUE,label_offset = 1) +} } diff --git a/man/assemble_genetic_features.Rd b/man/assemble_genetic_features.Rd new file mode 100644 index 0000000..ab1596b --- /dev/null +++ b/man/assemble_genetic_features.Rd @@ -0,0 +1,92 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{assemble_genetic_features} +\alias{assemble_genetic_features} +\title{Assemble genetic features for UMAP input} +\usage{ +assemble_genetic_features( + these_samples_metadata, + metadata_columns = c("bcl2_ba", "bcl6_ba", "myc_ba"), + genes, + synon_genes, + maf_with_synon, + hotspot_genes, + genome_build = "grch37", + sv_value = 3, + synon_value = 1, + coding_value = 2, + include_ashm = TRUE, + annotated_sv, + include_GAMBL_sv = TRUE, + review_hotspots = TRUE, + verbose = FALSE +) +} +\arguments{ +\item{these_samples_metadata}{Data frame with sample metadata, must include seq_type and sample_id.} + +\item{metadata_columns}{Columns in metadata to use for SV status (default: c("bcl2_ba","bcl6_ba","myc_ba")).} + +\item{genes}{Vector of gene symbols to include.} + +\item{synon_genes}{Vector of gene symbols for synonymous mutations.} + +\item{maf_with_synon}{MAF data frame including synonymous mutations.} + +\item{hotspot_genes}{Vector of hotspot genes.} + +\item{genome_build}{Genome build version (default: "grch37").} + +\item{sv_value}{Value to assign for SV presence (default: 3).} + +\item{synon_value}{Value to assign for synonymous mutations (default: 1).} + +\item{coding_value}{Value to assign for coding mutations (default: 2).} + +\item{include_ashm}{Logical; if TRUE, includes aSHM counts in the output matrix (default: TRUE).} + +\item{annotated_sv}{Data frame with annotated structural variants.} + +\item{include_GAMBL_sv}{Logical; if TRUE, includes GAMBL structural variants in the output matrix (default: TRUE).} + +\item{review_hotspots}{Logical; if TRUE, reviews hotspot genes (default: TRUE).} + +\item{verbose}{Defaults to FALSE} +} +\value{ +Matrix of assembled features for each sample. +} +\description{ +This function assembles a matrix of genetic features for each sample, including mutation status, +aSHM counts, and structural variant status for BCL2, BCL6, and MYC. It supports both genome and capture sequencing types. +} +\examples{ + +\dontrun{ +library(GAMBLR.predict) + +ordered_genes = filter(lymphoma_genes,DLBCL_Tier==1 | FL_Tier==1|BL_Tier==1) \%>\% pull(Gene) +synon_genes = unique(grch37_ashm_regions$gene) +synon_genes = synon_genes[synon_genes \%in\% ordered_genes] + +all_sv_anno = get_combined_sv(these_samples_metadata = dlbcl_meta) + +all_sv_anno = annotate_sv(all_sv_anno) + +all_full_status = assemble_genetic_features( + these_samples_metadata = dlbcl_meta, + metadata_columns = c("BCL2_SV","BCL6_SV"), + synon_genes = synon_genes, + synon_value = 1, + coding_value = 2, + genes = ordered_genes, + hotspot_genes = c("MYD88"), + maf_with_synon = all_maf_with_s, + include_ashm = F, + sv_value = 2, + verbose=F, + annotated_sv = all_sv_anno +) +} + +} diff --git a/man/basic_umap_scatterplot.Rd b/man/basic_umap_scatterplot.Rd new file mode 100644 index 0000000..7082a23 --- /dev/null +++ b/man/basic_umap_scatterplot.Rd @@ -0,0 +1,60 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{basic_umap_scatterplot} +\alias{basic_umap_scatterplot} +\title{Basic UMAP Scatterplot} +\usage{ +basic_umap_scatterplot( + optimized, + plot_samples = NULL, + colour_by = NULL, + truth_column = "lymphgen", + pred_column = "DLBCLone_ko", + other_label = "Other", + title = "UMAP based on selected features", + use_plotly = TRUE, + custom_colours = NULL +) +} +\arguments{ +\item{optimized}{Data frame containing at least V1, V2, sample_id, and grouping columns.} + +\item{plot_samples}{Optional character vector of sample_ids to label in the plot} + +\item{colour_by}{Column name to color points by. Defaults to \code{truth_column}.} + +\item{truth_column}{Name of the truth/ground-truth column (default: "lymphgen").} + +\item{pred_column}{Name of the predicted-class column (default: "DLBCLone_ko").} + +\item{other_label}{Label used for the outgroup/unclassified class (default: "Other").} + +\item{title}{Plot title.} + +\item{use_plotly}{Logical; if FALSE and \code{plot_samples} provided, draw static labels.} + +\item{custom_colours}{Optional named vector of colors for groups; falls back to \code{get_gambl_colours()}.} +} +\value{ +A ggplot object. +} +\description{ +Generates a simple UMAP scatterplot for visualizing sample clustering or separation. +} +\examples{ + +\dontrun{ +library(GAMBLR.predict) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) + +basic_umap_scatterplot( + make_umap$df, #the data frame containing V1 and V2 from UMAP + plot_samples = "some_sample_ID", + colour_by = "DLBCLone_ko") +} + +} diff --git a/man/construct_reduced_winning_version.Rd b/man/construct_reduced_winning_version.Rd index 6b66577..5e0772a 100644 --- a/man/construct_reduced_winning_version.Rd +++ b/man/construct_reduced_winning_version.Rd @@ -5,17 +5,19 @@ \title{Construct reduced 21-dimension feature vector for DLBCLass} \usage{ construct_reduced_winning_version( - mutations = "inst/extdata/DLBCL.699.fullGSM.Sep_23_2022.tsv", - fisher_test_result = "inst/extdata/fisher_exact_5x2.Sep_23_2022.combined.tsv", + mutations_file = "inst/extdata/DLBCL.699.fullGSM.Sep_23_2022.tsv", + fisher_test_result_file = "inst/extdata/fisher_exact_5x2.Sep_23_2022.combined.tsv", + include_cn = TRUE, + mutation_data, add_missing_features = FALSE ) } \arguments{ -\item{fisher_test_result}{Fisher's test result from DLBCLass github repository} - \item{add_missing_features}{Set to TRUE to fill in missing features with zeroes} \item{data}{Data frame of mutation status from DLBCLass supplement} + +\item{fisher_test_result}{Fisher's test result from DLBCLass github repository} } \value{ a list of data frames with the full mutation features, collapsed diff --git a/man/make_alluvial.Rd b/man/make_alluvial.Rd new file mode 100644 index 0000000..cac4e43 --- /dev/null +++ b/man/make_alluvial.Rd @@ -0,0 +1,101 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{make_alluvial} +\alias{make_alluvial} +\title{Create an Alluvial Plot Comparing Original and Predicted Classifications} +\usage{ +make_alluvial( + optimized, + count_excluded_as_other = FALSE, + title = "", + group_order = NULL, + add_accuracy_to_title = TRUE, + accuracy_per_group = TRUE, + accuracy_type = "sensitivity", + truth_name = "Lymphgen", + truth_column = "lymphgen", + pred_name = "DLBCLone", + pred_column = "predicted_label_optimized", + nudge = 0.03, + box_nudge = 0.15, + min_flip_n = 10, + add_percent = TRUE, + add_unclass_rate = TRUE, + denom = 20, + label_size = 3, + label_lock_y = TRUE, + label_line_flip_colour = TRUE, + label_box_flip_colour = TRUE, + concordant_label_relative_pos = 0.5, + verbose = FALSE +) +} +\arguments{ +\item{optimized}{List containing prediction results and metadata, typically output from a DLBCLone optimization function.} + +\item{count_excluded_as_other}{Logical; if TRUE, samples excluded due to missing features are counted as "Other" in the plot.} + +\item{title}{Plot title (default: empty string).} + +\item{group_order}{Character vector specifying the order of groups/classes for axes and coloring.} + +\item{add_accuracy_to_title}{Logical; if TRUE, adds accuracy/concordance rate to the plot title.} + +\item{accuracy_per_group}{Logical; if TRUE, computes and displays per-group accuracy.} + +\item{accuracy_type}{Type of accuracy to report (default: "sensitivity").} + +\item{truth_name}{Name for the original class column (default: "Lymphgen").} + +\item{truth_column}{Name for the original class column in the predictions data frame (default: "lymphgen").} + +\item{pred_name}{Name for the predicted class column (default: "DLBCLone").} + +\item{pred_column}{Name of the column in predictions to use for the predicted class (default: "predicted_label_optimized").} + +\item{nudge}{Amount to nudge labels horizontally (default: 0.03).} + +\item{box_nudge}{Amount to nudge label boxes (default: 0.15).} + +\item{min_flip_n}{Minimum number of samples for a flow to be labeled (default: 10).} + +\item{add_percent}{Logical; if TRUE, adds percent concordance to labels.} + +\item{add_unclass_rate}{Logical; if TRUE, adds unclassified rate annotation to the plot.} + +\item{denom}{Denominator for stratum/flow width scaling (default: 20).} + +\item{label_size}{Size of label text (default: 3).} + +\item{label_lock_y}{Logical; if TRUE, locks label movement to the x-axis only.} + +\item{label_line_flip_colour}{Logical; controls color assignment for label lines.} + +\item{label_box_flip_colour}{Logical; controls color assignment for label boxes.} + +\item{concordant_label_relative_pos}{Position for concordant labels: 0 (left), 0.5 (middle), or 1 (right).} + +\item{verbose}{Whether to print verbose outputs to console} +} +\value{ +A ggplot2 object representing the alluvial plot. +} +\description{ +This function generates a detailed alluvial plot to visualize the concordance and discordance between original (e.g., Lymphgen) and predicted (e.g., DLBCLone) class assignments for samples. +It supports annotation of concordance rates, per-group accuracy, unclassified rates, and flexible labeling and coloring options. +} +\details{ +\itemize{ +\item Visualizes flows between original and predicted classes, highlighting concordant and discordant assignments. +\item Annotates concordance rate, per-group accuracy, and unclassified rate as specified. +\item Supports flexible labeling, coloring, and axis ordering for publication-quality plots. +} +} +\examples{ +\dontrun{ +library(GAMBLR.predict) + +make_alluvial(optimize_params) +} + +} diff --git a/man/make_and_annotate_umap.Rd b/man/make_and_annotate_umap.Rd index 38a5d9e..a205a6c 100644 --- a/man/make_and_annotate_umap.Rd +++ b/man/make_and_annotate_umap.Rd @@ -8,6 +8,10 @@ make_and_annotate_umap( df, metadata, umap_out, + truth_column = "lymphgen", + core_features = NULL, + core_feature_multiplier = 1.5, + hidden_features = NULL, n_neighbors = 55, min_dist = 0, metric = "cosine", @@ -16,7 +20,13 @@ make_and_annotate_umap( ret_model = TRUE, na_vals = "drop", join_column = "sample_id", - seed = 42 + seed = 12345, + target_column, + target_metric = "euclidean", + target_weight = 0.5, + calc_dispersion = FALSE, + algorithm = "tumap", + make_plot = FALSE ) } \arguments{ @@ -29,6 +39,16 @@ matches the row names of df. This data frame will be joined to the UMAP output.} will use this model to project the data instead of re-running UMAP. This is useful for reproducibility and for using the same UMAP model on different datasets.} +\item{core_features}{A vector of column names in df that should be multiplied by +core_feature_multiplier. These features are considered "core" and will be weighted +more heavily in the UMAP embedding.} + +\item{core_feature_multiplier}{A numeric value that will be used to multiply the +core features specified in core_features. Default is 1.5.} + +\item{hidden_features}{A vector of column names in df that should be dropped from the +UMAP input. These features are considered "hidden" and will not be used in the UMAP embedding.} + \item{n_neighbors}{Passed to UMAP2. The number of neighbors to consider when calculating the UMAP embedding.} \item{min_dist}{Passed to UMAP2. The minimum distance between points in the UMAP embedding.} @@ -48,21 +68,34 @@ all NA to zero and leaves the column intact.} \item{join_column}{The column name in the metadata data frame that contains the sample IDs (default sample_id).} \item{seed}{Passed to UMAP2. The random seed for reproducibility.} + +\item{target_column}{The column in the metadata that contains the target variable for classification.} + +\item{target_metric}{The distance metric to use for calculating distances to the target variable.} + +\item{target_weight}{The weight to assign to the target variable in the UMAP embedding. Default is 0.5.} + +\item{calc_dispersion}{If TRUE, calculates the dispersion of the UMAP embedding.} + +\item{algorithm}{The UMAP algorithm to use. Default is "tumap", which uses the TUMAP implementation.} + +\item{make_plot}{If TRUE, creates a plot of the UMAP embedding.} } \description{ Run UMAP and attach result to metadata } \examples{ -umap_outs = make_and_annotate_umap(df=gambl_train_lymphgen, - min_dist = 0, - n_neighbors = 55, - init="spectral", - n_epochs = 1500, - #seed=best_params$seed, - metadata=gambl_train_meta_dlbclass, - ret_model=T, - metric="cosine") +library(GAMBLR.predict) +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) } diff --git a/man/make_neighborhood_plot.Rd b/man/make_neighborhood_plot.Rd new file mode 100644 index 0000000..37cabb7 --- /dev/null +++ b/man/make_neighborhood_plot.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{make_neighborhood_plot} +\alias{make_neighborhood_plot} +\title{Make Neighborhood Plot} +\usage{ +make_neighborhood_plot( + single_sample_prediction_output, + training_predictions, + this_sample_id, + prediction_in_title = TRUE, + add_circle = TRUE, + label_column = "DLBCLone_io", + point_size = 0.5, + point_alpha = 0.9, + line_alpha = 0.9 +) +} +\arguments{ +\item{single_sample_prediction_output}{A list containing prediction results and annotation data frames. +Must include elements \code{prediction} (data frame with prediction results) and \code{anno_df} (data frame with UMAP coordinates and annotations).} + +\item{training_predictions}{The equivalent data frame of prediction results for all training samples (e.g. optimized_model$df)} + +\item{this_sample_id}{Character. The sample ID for which the neighborhood plot will be generated.} + +\item{prediction_in_title}{Logical. If \code{TRUE}, includes the predicted label in the plot title.} + +\item{add_circle}{Plot will include a circle surrounding the set of neighbors. Set to FALSE to disable.} + +\item{label_column}{Does nothing, i.e. this is not currently working.} +} +\value{ +A \code{ggplot2} object representing the UMAP plot with the selected sample and its neighbors highlighted. +} +\description{ +Generates a UMAP plot highlighting the neighborhood of a given sample, showing its nearest neighbors and their group assignments. +} +\details{ +The function extracts the nearest neighbors of the specified sample, draws segments connecting the sample to its neighbors, and colors points by group (e.g., lymphgen subtype). The plot title can optionally include the predicted label. +} +\examples{ + +# Assuming 'optimization_result' is the output of DLBCLone_optimize_params +# and 'output' is the result of DLBCLone_predict_single_sample +# on sample_id "SAMPLE123": +\dontrun{ +library(GAMBLR.predict) + +predict_single <- predict_single_sample_DLBCLone( + test_df = optimize_params$features[1,], + train_metadata = dlbcl_meta, + optimized_model = optimize_params +) + +make_neighborhood_plot( + single_sample_prediction_output = predict_single, + training_predictions = optimize_params$df, + this_sample_id = "SAMPLE123", + prediction_in_title = TRUE, + add_circle = TRUE +) +} + +} diff --git a/man/make_umap_scatterplot.Rd b/man/make_umap_scatterplot.Rd new file mode 100644 index 0000000..7c29ffb --- /dev/null +++ b/man/make_umap_scatterplot.Rd @@ -0,0 +1,56 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{make_umap_scatterplot} +\alias{make_umap_scatterplot} +\title{Make UMAP scatterplot} +\usage{ +make_umap_scatterplot( + df, + truth_column = "lymphgen", + drop_composite = TRUE, + drop_other = FALSE, + high_confidence = FALSE, + custom_colours, + add_labels = FALSE, + title = NULL +) +} +\arguments{ +\item{df}{Data frame containing the UMAP coordinates and annotations.} + +\item{truth_column}{Name of the column containing the labels to be plotted (default: "lymphgen").} + +\item{drop_composite}{If TRUE: removes composite labels from the lymphgen column.} + +\item{drop_other}{If TRUE: removes "Other" and "NOS" labels from the lymphgen column.} + +\item{high_confidence}{If TRUE: filters the data to include only samples with confidence > 0.7.} + +\item{custom_colours}{Custom color palette for the plot. If not provided, uses default GAMBL colors.} + +\item{add_labels}{If TRUE: adds labels to the points based on the median coordinates of each group.} + +\item{title}{Title for the plot. Default: NULL.} +} +\value{ +A ggplot object representing the UMAP scatterplot with marginal histograms. +} +\description{ +Make UMAP scatterplot +} +\examples{ +\dontrun{ +library(GAMBLR.predict) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) + +make_umap_scatterplot( + make_umap$df, + drop_other = F +) +} + +} diff --git a/man/nearest_neighbor_heatmap.Rd b/man/nearest_neighbor_heatmap.Rd new file mode 100644 index 0000000..17e4858 --- /dev/null +++ b/man/nearest_neighbor_heatmap.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{nearest_neighbor_heatmap} +\alias{nearest_neighbor_heatmap} +\title{Heatmap visualization of mutations in nearest neighbors for a sample} +\usage{ +nearest_neighbor_heatmap( + this_sample_id, + DLBCLone_model, + truth_column = "lymphgen", + metadata_cols = NULL, + clustering_distance = "binary", + font_size = 14 +) +} +\arguments{ +\item{this_sample_id}{Character. The sample ID for which to plot the nearest neighbor heatmap.} + +\item{DLBCLone_model}{A DLBCLone model object, which can be the output of \code{DLBCLone_optimize_params}, \code{DLBCLone_KNN}, or \code{predict_single_sample_DLBCLone}.} + +\item{truth_column}{The column name in the predictions data frame that contains the ground truth labels (default: "lymphgen").} + +\item{metadata_cols}{Optional character vector of additional metadata columns to include in the heatmap annotations.} + +\item{clustering_distance}{The distance metric to use for clustering rows. Default is "binary".} + +\item{font_size}{Font size for heatmap text (default: 14).} +} +\value{ +A ComplexHeatmap object showing the feature matrix for the nearest neighbors of the sample. +} +\description{ +Generates a heatmap of feature values for the nearest neighbors of a specified sample, +based on a DLBCLone model object. This visualization helps to inspect the feature profiles +of samples most similar to the query sample. +} +\details{ +\itemize{ +\item For models of type \code{DLBCLone_optimize_params}, uses the \code{neighbors} and \code{lyseq_status} fields. +\item For models of type \code{DLBCLone_KNN}, uses the \code{unlabeled_neighbors} field. +\item The function extracts the feature matrix rows corresponding to the nearest neighbors of the specified sample, +and plots a heatmap of features with nonzero values. +} +} +\examples{ +# Assuming 'model' is a DLBCLone model and 'sample_id' is a valid sample: +\dontrun{ +library(GAMBLR.predict) + +nearest_neighbor_heatmap( + this_sample_id = "CCS_0680_lst", + DLBCLone_model = predict_single, + font_size = 10 +) +} + +} diff --git a/man/optimize_outgroup.Rd b/man/optimize_outgroup.Rd new file mode 100644 index 0000000..d7dd9bc --- /dev/null +++ b/man/optimize_outgroup.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{optimize_outgroup} +\alias{optimize_outgroup} +\title{Optimize the threshold for classifying samples as "Other"} +\usage{ +optimize_outgroup( + predicted_labels, + true_labels, + other_score, + all_classes = c("MCD", "EZB", "BN2", "N1", "ST2", "Other"), + maximize = "balanced_accuracy", + exclude_other_for_accuracy = FALSE, + cap_classification_rate = 1, + verbose = FALSE, + other_class = "Other" +) +} +\arguments{ +\item{predicted_labels}{Vector of predicted labels for the samples} + +\item{true_labels}{Vector of true labels for the samples} + +\item{other_score}{Vector of scores for the "Other" class for each sample ()} + +\item{all_classes}{Vector of classes to use for training and testing. +Default: c("MCD","EZB","BN2","N1","ST2","Other")} + +\item{maximize}{Metric to use for optimization. +Either "accuracy" (actual accuracy across all samples) or "balanced_accuracy" +(the mean of the balanced accuracy values across all classes). +Default: "balanced_accuracy"} + +\item{exclude_other_for_accuracy}{Set to TRUE to exclude the +"Other" class from the 'lymphgen' column when calculating accuracy metrics +(passed to DLBCLone_optimize_params). Default: FALSE} + +\item{verbose}{Set to TRUE to print additional information during the optimization process.} +} +\value{ +a list of data frames with the predictions and the UMAP input +} +\description{ +Performs a post-hoc evaluation of the classification of a sample as one of +the main classes vs the outgroup/unclassified label "Other" and returns the +optimal threshold for classifying a sample as "Other" based on the ground +truth provided in the true_labels vector. It evaluates the performance +of the classifier using a range of thresholds and returns the best threshold +based on the specified metric (balanced accuracy or accuracy). +This function is not generally meant to be called directly but rather is +a helper function used by DLBCLone_optimize_params. +} diff --git a/man/optimize_purity.Rd b/man/optimize_purity.Rd new file mode 100644 index 0000000..f453d39 --- /dev/null +++ b/man/optimize_purity.Rd @@ -0,0 +1,60 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{optimize_purity} +\alias{optimize_purity} +\title{Optimize Purity Threshold for Classification Assignment} +\usage{ +optimize_purity( + optimized_model_object, + vote_df, + mode, + optimize_by = "balanced_accuracy", + truth_column, + all_classes = c("MCD", "EZB", "BN2", "N1", "ST2", "Other"), + k, + cap_classification_rate = 1, + exclude_other_for_accuracy = FALSE, + other_class = "Other" +) +} +\arguments{ +\item{optimized_model_object}{An object containing the optimized model parameters and predictions.} + +\item{vote_df}{Data frame containing the kNN vote results, including columns for sample IDs, predicted labels, and scores.} + +\item{mode}{The mode of operation, which determines the output column names (e.g., "DLBCLone_w").} + +\item{optimize_by}{The metric to optimize by, either "balanced_accuracy"} + +\item{truth_column}{Name of the column in \code{processed_votes} containing the true class labels.} + +\item{all_classes}{Vector of all possible class labels to consider (default: c("MCD", "EZB", "BN2", "N1", "ST2", "Other")).} + +\item{k}{The number of neighbors used in the kNN classification.} + +\item{exclude_other_for_accuracy}{Logical indicating whether to exclude the "Other" class from accuracy calculations (default: FALSE).} +} +\value{ +A list with two elements: \code{best_accuracy} (numeric, the highest accuracy achieved) and \code{best_purity_threshold} (numeric, the threshold at which this accuracy was achieved). +} +\description{ +This function searches for the optimal purity threshold to assign samples to their predicted class or to "Other" based on the score ratio in processed kNN vote results. +It iteratively tests a range of purity thresholds, updating the predicted class if the score ratio meets or exceeds the threshold, and computes the accuracy for each threshold. +The function returns the best accuracy achieved and the corresponding purity threshold. +} +\details{ +\itemize{ +\item For each threshold in the range 0.1 to 0.95 (step 0.05), the function updates the prediction column to assign the class from \code{by_score_opt} if the score ratio meets the threshold, otherwise assigns "Other". +\item Accuracy is computed as the proportion of correct assignments (diagonal of the confusion matrix). +\item The function is intended for use in optimizing classification purity in kNN-based workflows, especially when distinguishing between confident class assignments and ambiguous ("Other") cases. +} +} +\examples{ + +\dontrun{ +library(GAMBLR.predict) + +result <- optimize_purity(processed_votes, prediction_column = "pred_label", truth_column = "true_label") +} + +} diff --git a/man/predict_single_sample.Rd b/man/predict_single_sample.Rd deleted file mode 100644 index e429880..0000000 --- a/man/predict_single_sample.Rd +++ /dev/null @@ -1,37 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/umap.R -\name{predict_single_sample} -\alias{predict_single_sample} -\title{Predict class for a single sample without using umap_transform} -\usage{ -predict_single_sample( - test_df, - train_df, - train_metadata, - best_params, - truth_classes = c("EZB", "MCD", "ST2", "N1", "BN2"), - drop_unlabeled_from_training = TRUE, - make_plot = TRUE -) -} -\arguments{ -\item{test_df}{Data frame containing the mutation status of the test sample} - -\item{train_df}{Data frame containing the mutation status of the training samples} - -\item{train_metadata}{Metadata for training samples with truth labels in lymphgen column} - -\item{best_params}{Data frame from DLBCLone_optimize_params with the best parameters} - -\item{truth_classes}{Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2")} - -\item{drop_unlabeled_from_training}{Set to TRUE to drop unlabeled samples from the training data} - -\item{make_plot}{Set to TRUE to plot the UMAP projection and predictions} -} -\value{ -a list of data frames with the predictions and the UMAP input -} -\description{ -Predict class for a single sample without using umap_transform -} diff --git a/man/predict_single_sample_DLBCLone.Rd b/man/predict_single_sample_DLBCLone.Rd new file mode 100644 index 0000000..b0fa4b4 --- /dev/null +++ b/man/predict_single_sample_DLBCLone.Rd @@ -0,0 +1,66 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{predict_single_sample_DLBCLone} +\alias{predict_single_sample_DLBCLone} +\title{Predict class for a single sample without using umap_transform and plot result of classification} +\usage{ +predict_single_sample_DLBCLone( + test_df, + train_metadata, + optimized_model = NULL, + truth_classes = c("EZB", "MCD", "ST2", "N1", "BN2", "Other"), + seed = 12345, + max_neighbors = 500 +) +} +\arguments{ +\item{test_df}{Data frame containing the mutation status of the test sample} + +\item{train_metadata}{Metadata for training samples with truth labels in lymphgen column} + +\item{optimized_model}{list of parameters from DLBCLone_optimize_params, neccessary UMAP output from a +previous, Data frame with the best parameters. useful for reproducibility.} + +\item{truth_classes}{Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2")} + +\item{seed}{Random seed for reproducibility} + +\item{max_neighbors}{Maximum number of neighbors to consider for each sample. Default 500.} +} +\value{ +a list of data frames with the predictions, the UMAP input, the model, and a ggplot object +} +\description{ +Predict class for a single sample without using umap_transform and plot result of classification +} +\examples{ + +library(GAMBLR.predict) + +dlbcl_meta = readr::read_tsv(system.file("extdata/dlbcl_meta_with_dlbclass.tsv",package = "GAMBLR.predict")) + +all_full_status = readr::read_tsv(system.file("extdata/all_full_status.tsv",package = "GAMBLR.predict")) \%>\% + tibble::column_to_rownames("sample_id") + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) + +optimize_params <- DLBCLone_optimize_params( + all_full_status, + dlbcl_meta, + umap_out = make_umap, + truth_classes = c("MCD","EZB","BN2","N1","ST2","Other"), + optimize_for_other = T, + min_k=5, + max_k=23 +) + +predict_single <- predict_single_sample_DLBCLone( + test_df = optimize_params$features[1,], + train_metadata = dlbcl_meta, + optimized_model = optimize_params +) + +} diff --git a/man/process_votes.Rd b/man/process_votes.Rd new file mode 100644 index 0000000..7a91863 --- /dev/null +++ b/man/process_votes.Rd @@ -0,0 +1,59 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{process_votes} +\alias{process_votes} +\title{Process KNN Vote Strings and Scores for Classification} +\usage{ +process_votes( + df, + raw_col = "label", + group_labels = c("EZB", "MCD", "ST2", "BN2", "N1", "Other"), + vote_labels_col = "vote_labels", + k, + other_vote_multiplier = 2, + score_purity_requirement = 1, + weighted_votes_col = "weighted_votes", + other_class = "Other" +) +} +\arguments{ +\item{df}{Data frame containing kNN results, including columns with neighbor labels and weighted votes.} + +\item{raw_col}{Name of the column containing the comma-separated neighbor labels (default: "label").} + +\item{group_labels}{Character vector of all possible class labels to consider (default: c("EZB", "MCD", "ST2", "BN2", "N1", "Other")).} + +\item{vote_labels_col}{Name of the column containing the comma-separated neighbor labels for weighted votes (default: "vote_labels").} + +\item{k}{Number of neighbors used in kNN (required).} + +\item{other_vote_multiplier}{Multiplier for the "Other" class when determining if a sample should be reclassified as "Other" (default: 2).} + +\item{score_purity_requirement}{Minimum ratio of top group score to "Other" score to assign a sample to the top group (default: 1).} + +\item{weighted_votes_col}{Name of the column containing the comma-separated weighted votes (default: "weighted_votes").} +} +\value{ +Data frame with additional columns for per-class neighbor counts, scores, top group assignments, and summary statistics for each sample. +} +\description{ +This function processes the raw neighbor label strings and weighted vote scores from k-nearest neighbor (KNN) classification results. +It computes per-class neighbor counts, weighted scores, and identifies the top group by count and score for each sample. +The function also supports custom logic for handling the "Other" class, including vote multipliers and purity requirements. +} +\details{ +\itemize{ +\item Computes the number of neighbors for each class and the sum of weighted votes per class. +\item Identifies the top group by count and by weighted score, and applies custom logic for the "Other" class if present. +\item Adds columns for counts, scores, top group, top group score, score ratios, and optimized group assignments. +\item Designed for downstream use in DLBCLone and similar kNN-based classification workflows. +} +} +\examples{ +/dontrun{ +library(GAMBLR.predict) + +result <- process_votes(knn_output_df, k = 7) +} + +} diff --git a/man/report_accuracy.Rd b/man/report_accuracy.Rd new file mode 100644 index 0000000..f9bdac1 --- /dev/null +++ b/man/report_accuracy.Rd @@ -0,0 +1,56 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{report_accuracy} +\alias{report_accuracy} +\title{Report Classification Accuracy and Per-Class Metrics} +\usage{ +report_accuracy( + predictions, + truth = "lymphgen", + pred = "DLBCLone_io", + per_group = FALSE, + metric = "accuracy", + verbose = FALSE +) +} +\arguments{ +\item{predictions}{Data frame containing predicted and true class labels.} + +\item{truth}{Name of the column with true class labels (default: "lymphgen").} + +\item{pred}{Name of the column with predicted class labels (default: "predicted_label").} + +\item{per_group}{Logical; if TRUE, computes per-group accuracy metrics.} + +\item{metric}{Character; type of accuracy to report ("accuracy" supported).} + +\item{verbose}{Whether to print verbose outputs to console} +} +\value{ +A list with: +\item{no_other}{Accuracy excluding samples assigned to "Other"} +\item{per_class}{Balanced accuracy per class} +\item{per_class_sensitivity}{Sensitivity per class} +\item{overall}{Overall accuracy including all samples} +} +\description{ +Computes overall accuracy, balanced accuracy, and sensitivity for predicted vs. true class labels. +Optionally excludes samples assigned to the "Other" class from accuracy calculations. +} +\details{ +\itemize{ +\item Uses confusion matrices to compute accuracy metrics. +\item Excludes "Other" class for no_other accuracy. +\item Returns per-class metrics for further analysis. +} +} +\examples{ +\dontrun{ +library(GAMBLR.predict) + +result <- report_accuracy(predictions_df) +result$overall +result$per_class +} + +} diff --git a/man/stacked_bar_plot.Rd b/man/stacked_bar_plot.Rd new file mode 100644 index 0000000..3c86438 --- /dev/null +++ b/man/stacked_bar_plot.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotting.R +\name{stacked_bar_plot} +\alias{stacked_bar_plot} +\title{Create a stacked bar plot of top features per subtype} +\usage{ +stacked_bar_plot( + DLBCLone_model, + truth_column = "lymphgen", + truth_classes = c("BN2", "EZB", "MCD", "ST2"), + method = "common", + num_feats = 10, + title = NULL +) +} +\arguments{ +\item{DLBCLone_model}{A DLBCLone model object, which can be the output of \code{DLBCLone_optimize_params}, or \code{DLBCLone_KNN}} + +\item{truth_column}{Name of the column containing the true class labels (default: "lymphgen").} + +\item{truth_classes}{Vector of class labels to consider (default: c("BN2","EZB","MCD","ST2")).} + +\item{method}{Method to determine top features: "common" for most common features, "chi_square" for subtype vs rest significance (default : "common").} + +\item{num_feats}{Number of top features to display per subtype (default: 10).} + +\item{title}{Title for the plot (default: NULL).} +} +\value{ +A ggplot2 object representing the stacked bar plot. +} +\description{ +This function generates a stacked bar plot showing the top features (genes) for each subtype based on their prevalence in the dataset. +} +\examples{ +\dontrun{ +library(GAMBLR.predict) + +stacked_bar_plot( + optimize_params, + method = "chi_square", + num_feats = 10, + title = "LymphGen" +) +} + +} diff --git a/man/summarize_all_ssm_status.Rd b/man/summarize_all_ssm_status.Rd new file mode 100644 index 0000000..98f48a0 --- /dev/null +++ b/man/summarize_all_ssm_status.Rd @@ -0,0 +1,86 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/umap.R +\name{summarize_all_ssm_status} +\alias{summarize_all_ssm_status} +\title{Summarize SSM (Somatic Single Nucleotide Mutation) Status Across Samples} +\usage{ +summarize_all_ssm_status( + maf_df, + these_samples_metadata, + genes_of_interest, + synon_genes, + silent_maf_df, + separate_by_class_genes = NULL, + count_hits = FALSE +) +} +\arguments{ +\item{maf_df}{A data frame containing mutation annotation format (MAF) data, +with at least the following columns: +\code{Hugo_Symbol}, \code{Variant_Classification}, and \code{Tumor_Sample_Barcode}.} + +\item{these_samples_metadata}{A data frame containing metadata for the samples, +with at least a \code{sample_id} column. +Any sample that does not have a matching sample_id in these_samples_metadata will be dropped.} + +\item{genes_of_interest}{A character vector of gene symbols to include +in the summary. If missing, defaults to all Tier 1 B-cell lymphoma genes.} + +\item{synon_genes}{(Optional) A character vector of gene symbols for which +synonymous mutations should be included.} + +\item{silent_maf_df}{(Optional) A separate data frame containing silent mutation data if +the user doesn't want to pull silent mutation status from \code{maf_df}. +This argument is useful when you want to combine mutations from +the output of get_coding_ssm and get_ssm_by_region or get_ssm_by_gene} + +\item{separate_by_class_genes}{(Optional) A character vector of +gene symbols for which mutations should be separated by class +(e.g., "Nonsense_Mutation", "Missense_Mutation").} + +\item{count_hits}{Logical; if \code{TRUE}, counts the number of mutations +per gene per sample. If \code{FALSE} (default), only presence/absence is recorded.} +} +\value{ +A wide-format data frame (matrix) with samples as rows and +mutation types as columns. Each cell contains either the count of +mutations (if \code{count_hits = TRUE}) or a binary indicator (0/1) +for mutation presence. +} +\description{ +This function summarizes the mutation status for a set of genes +across multiple samples, separating mutations by class for genes specified +by \code{separate_by_class_genes} and counting the number of hits +in each sample per mutation category. +} +\details{ +\itemize{ +\item Mutations are grouped and optionally separated by mutation +class for genes specified in \code{separate_by_class_genes} +\item Synonymous mutations can be counted as another separate feature for genes specified by \code{synon_genes} genes. +\item The function simplifies mutation annotations and pivots the data to a wide format suitable for downstream analysis. +} +} +\examples{ +# A basic example, using only the output of get_all_coding_ssm +# Since the only non-coding class this function handles is Silent, +# we will be missing most non-coding types such as Intron, UTR, Flank + +\dontrun{ +library(GAMBLR.predict) + +sample_metadata = get_sample_metadata() \%>\% filter(seq_type!= "mrna") + +maf_data = get_all_coding_ssm(sample_metadata) + +mutation_matrix <- summarize_all_ssm_status( + maf_df = maf_data, + these_samples_metadata = sample_metadata, + genes_of_interest = c("TP53", "SGK1", "BCL2"), + synon_genes = c("BCL2"), + separate_by_class_genes = c("TP53","SGK1"), + count_hits = FALSE +) +} + +} diff --git a/man/weighted_knn_predict_with_conf.Rd b/man/weighted_knn_predict_with_conf.Rd index 21821ef..4c1afc5 100644 --- a/man/weighted_knn_predict_with_conf.Rd +++ b/man/weighted_knn_predict_with_conf.Rd @@ -9,13 +9,15 @@ weighted_knn_predict_with_conf( train_labels, test_coords, k, - epsilon = 1e-05, + epsilon = 0.1, conf_threshold = NULL, - na_label = "Other", + other_class = "Other", verbose = FALSE, use_weights = TRUE, - ignore_top = FALSE, - track_neighbors = FALSE + ignore_self = TRUE, + track_neighbors = TRUE, + separate_other = TRUE, + max_neighbors = 500 ) } \arguments{ @@ -36,22 +38,37 @@ when use_weights is TRUE. Default: 1e-5} \item{conf_threshold}{Minimum confidence for classifying a sample based on it's neighbors. Below this value the sample will be assigned na_label instead} -\item{na_label}{Class to assign all samples that are not confidently assigned. -Default: Other} - \item{verbose}{Whether to print verbose outputs to console} \item{use_weights}{Set to FALSE for all neigbors to have equal weight when calculating the confidence} -\item{ignore_top}{Set to TRUE to avoid considering a nearest neighbor with -distance = 0. This is usually only relevant when re-classifying labeled -samples to estimate overall accuracy} - \item{track_neighbors}{Set to TRUE to include details for the nearest neighbors of each sample} + +\item{separate_other}{Set to TRUE to treat the "Other" class separately +when calculating the confidence.} + +\item{max_neighbors}{Maximum number of neighbors to consider for each sample. Default 500.} + +\item{na_label}{Class to assign all samples that are not confidently assigned. +Default: Other} } \value{ -data frame with labels and confidence values for rows in test_coords +Data frame with rows = test samples and columns: +\item{predicted_label}{the predicted class} +\item{confidence}{predicted class weight / total weight} +If \code{track_neighbors = TRUE}, additional columns: +\item{other_score}{relative weight of outgroup vs predicted class} +\item{neighbor_id}{comma-separated neighbor sample IDs} +\item{neighbor}{comma-separated neighbor indices (in train order)} +\item{distance}{comma-separated neighbor distances} +\item{label}{comma-separated neighbor labels (in-group only if separate_other=TRUE)} +\item{vote_labels}{comma-separated unique labels contributing to weights} +\item{weighted_votes}{comma-separated weights per \code{vote_labels}} +\item{neighbors_other}{count of outgroup neighbors closer than the farthest in-group neighbor} +\item{other_weighted_votes}{sum of outgroup weights closer than the farthest in-group neighbor} +\item{total_w}{sum of weights for in-group neighbors used} +\item{pred_w}{weight supporting the predicted class} } \description{ Weighted k-nearest neighbor with confidence estimate