From ad1782da06cf4073dab2edb51e1e2728a1ffe83c Mon Sep 17 00:00:00 2001 From: Atlantis-Real Date: Tue, 16 Sep 2025 23:28:17 -0700 Subject: [PATCH 1/5] adding stacked_barplot to plotting.R --- NAMESPACE | 2 + R/plotting.R | 173 ++++++++++++++++++++++++++++++++++++++++ man/stacked_bar_plot.Rd | 47 +++++++++++ 3 files changed, 222 insertions(+) create mode 100644 man/stacked_bar_plot.Rd diff --git a/NAMESPACE b/NAMESPACE index 3739b01..c1f8c81 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -31,6 +31,7 @@ export(optimize_purity) 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) @@ -50,6 +51,7 @@ import(plotly) import(purrr) import(randomForest, except = c("combine")) import(readr) +import(rlang) import(shiny) import(shinybusy) import(shinyjs) diff --git a/R/plotting.R b/R/plotting.R index 8c7b6cf..c52a0fc 100644 --- a/R/plotting.R +++ b/R/plotting.R @@ -1642,3 +1642,176 @@ if(add_percent){ segment.colour = "none", fill = guide_legend(title = "Class")) } + +#' 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/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" +) +} + +} From 69949aa7a4cafd361e4214c62027f288aea25af8 Mon Sep 17 00:00:00 2001 From: Atlantis-Real Date: Tue, 16 Sep 2025 23:44:06 -0700 Subject: [PATCH 2/5] documentation and documentation examples added to plotting.R as well as hashed out suggestions for make_neighborhood_plot --- NAMESPACE | 14 +-- R/plotting.R | 145 ++++++++++++++++++++------------ man/DLBCLone_summarize_model.Rd | 4 +- man/basic_umap_scatterplot.Rd | 18 ++-- man/make_alluvial.Rd | 7 +- man/make_neighborhood_plot.Rd | 17 +++- man/make_umap_scatterplot.Rd | 33 ++++++-- man/nearest_neighbor_heatmap.Rd | 11 ++- man/report_accuracy.Rd | 4 + 9 files changed, 170 insertions(+), 83 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index c1f8c81..549bb97 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -24,7 +24,6 @@ export(make_and_annotate_umap) export(make_neighborhood_plot) export(make_umap_scatterplot) export(massage_matrix_for_clustering) -export(nearest_neighbor_heatmap) export(og_nearest_neighbor_heatmap) export(optimize_outgroup) export(optimize_purity) @@ -41,11 +40,13 @@ 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(plotly) import(purrr) @@ -61,14 +62,5 @@ import(tidyselect) import(uwot) importFrom(circlize,colorRamp2) importFrom(dplyr,filter) -importFrom(dplyr,left_join) -importFrom(dplyr,mutate) -importFrom(dplyr,pull) -importFrom(dplyr,select) importFrom(grid,gpar) importFrom(magrittr,"%>%") -importFrom(rlang,sym) -importFrom(tibble,column_to_rownames) -importFrom(tibble,rownames_to_column) -importFrom(tidyr,pivot_longer) -importFrom(tidyr,separate) diff --git a/R/plotting.R b/R/plotting.R index c52a0fc..4ae2b84 100644 --- a/R/plotting.R +++ b/R/plotting.R @@ -14,16 +14,21 @@ #' @param font_size Numeric. Font size for labels (default 14). #' #' @return A ComplexHeatmap object (drawn). -#' @importFrom dplyr filter select left_join mutate pull -#' @importFrom tidyr separate pivot_longer -#' @importFrom tibble rownames_to_column column_to_rownames -#' @import ComplexHeatmap -#' @importFrom grid gpar -#' @importFrom circlize colorRamp2 +#' +#' @import dplyr tidyr ComplexHeatmap grid circlize tibble +#' #' @examples -#' # Assuming 'predicted_out' is a the output of DLBCLone_KNN_predict +#' # Assuming 'model' is a DLBCLone model and 'sample_id' is a valid sample: +#' \dontrun{ +#' library(GAMBLR.predict) #' -#' @export +#' 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, @@ -300,10 +305,6 @@ nearest_neighbor_heatmap <- function( } - - - - ComplexHeatmap::draw( ht, heatmap_legend_side = "bottom", @@ -311,8 +312,6 @@ nearest_neighbor_heatmap <- function( ) } - - #' Heatmap visualization of mutations in nearest neighbors for a sample #' #' Generates a heatmap of feature values for the nearest neighbors of a specified sample, @@ -343,8 +342,6 @@ nearest_neighbor_heatmap <- function( #' nearest_neighbor_heatmap("SomeSample_ID", predicted_out) #'} #' - - depr_nearest_neighbor_heatmap <- function( this_sample_id, DLBCLone_model, @@ -691,18 +688,25 @@ og_nearest_neighbor_heatmap <- function(this_sample_id, #' @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{ -#' my_umap = make_and_annotate_umap(my_data, my_metadata) +#' library(GAMBLR.predict) +#' +#' make_umap <- make_and_annotate_umap( +#' df=all_full_status, +#' metadata=dlbcl_meta +#' ) #' -#' basic_umap_scatterplot(my_umap$df, #the data frame containing V1 and V2 from UMAP -#' plot_samples = "some_sample_ID", -#' colour_by = "DLBCLone_ko") +#' 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, @@ -781,7 +785,6 @@ message("colour_by: ", colour_by) return(p) } - #' Summarize and Export DLBCLone Model Results #' #' Generates and saves a set of summary plots and tables for @@ -803,15 +806,14 @@ message("colour_by: ", colour_by) #' #' @return No return value. Side effect: writes multiple PDF files to disk. #' -#' @import ggalluvial -#' +#' @import dplyr ggplot2 ComplexHeatmap rlang #' @export #' #' @examples #' \dontrun{ -#' DLBCLone_summarize_model(optimized_model = all_features_optimized, -#' base_name="model_summaries/Lymphgen_all_features") -#' } +#' DLBCLone_summarize_model("Full_geneset_unweighted", optimized_model) +#'} +#' DLBCLone_summarize_model = function(base_name, optimized_model){ base_dir = here::here() @@ -875,7 +877,6 @@ DLBCLone_summarize_model = function(base_name, } - #' @title Make Neighborhood Plot #' @description #' Generates a UMAP plot highlighting the neighborhood of a given sample, showing its nearest neighbors and their group assignments. @@ -893,19 +894,32 @@ DLBCLone_summarize_model = function(base_name, #' @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 -#' @import ggplot2 -#' @importFrom rlang sym -#' +#' @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{ -#' make_neighborhood_plot(output, optimization_result$df, "SAMPLE123") +#' 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, @@ -966,6 +980,7 @@ xmin = min(training_predictions$V1, na.rm = TRUE) } 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 = abs(V1-my_x)+abs(V2-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), @@ -987,26 +1002,45 @@ xmin = min(training_predictions$V1, na.rm = TRUE) 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) } +# 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 -#' @param drop_composite -#' @param colour_by -#' @param drop_other -#' @param high_confidence -#' @param custom_colours -#' @param add_labels +#' @param df Data frame containing the UMAP coordinates and annotations. +#' @param drop_composite If TRUE: removes composite labels from the lymphgen column. +#' @param colour_by Column name to color points by. Default: "lymphgen". +#' @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. #' -#' @returns +#' @returns A ggplot object representing the UMAP scatterplot with marginal histograms. #' -#' @import ggside +#' @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, drop_composite = TRUE, colour_by="lymphgen", @@ -1081,8 +1115,6 @@ make_umap_scatterplot = function(df, return(p) } - - #' Calculate Classification Accuracy and Per-Class Metrics based on Predictions #' #' Computes overall accuracy, balanced accuracy, and sensitivity for predicted vs. true class labels. @@ -1105,12 +1137,18 @@ make_umap_scatterplot = function(df, #' - 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 +#' } #' -#' @export report_accuracy <- function(predictions, truth = "lymphgen", pred = "DLBCLone_io", @@ -1211,7 +1249,6 @@ report_accuracy <- function(predictions, )) } - #' Create an Alluvial Plot Comparing Original and Predicted Classifications #' #' This function generates a detailed alluvial plot to visualize the concordance @@ -1253,12 +1290,16 @@ report_accuracy <- function(predictions, #' - Annotates concordance rate, per-group accuracy, and unclassified rate as specified. #' - Supports flexible labeling, coloring, and axis ordering for publication-quality plots. #' -#' @examples -#' # Example usage: -#' # make_alluvial(optimized_result) -#' -#' @import ggrepel +#' @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, diff --git a/man/DLBCLone_summarize_model.Rd b/man/DLBCLone_summarize_model.Rd index f57b91a..d8ac8d9 100644 --- a/man/DLBCLone_summarize_model.Rd +++ b/man/DLBCLone_summarize_model.Rd @@ -33,7 +33,7 @@ base name. } \examples{ \dontrun{ -DLBCLone_summarize_model(optimized_model = all_features_optimized, - base_name="model_summaries/Lymphgen_all_features") +DLBCLone_summarize_model("Full_geneset_unweighted", optimized_model) } + } diff --git a/man/basic_umap_scatterplot.Rd b/man/basic_umap_scatterplot.Rd index 86a2d3e..e0930de 100644 --- a/man/basic_umap_scatterplot.Rd +++ b/man/basic_umap_scatterplot.Rd @@ -35,19 +35,23 @@ basic_umap_scatterplot( \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{ -my_umap = make_and_annotate_umap(my_data, my_metadata) +library(GAMBLR.predict) + +make_umap <- make_and_annotate_umap( + df=all_full_status, + metadata=dlbcl_meta +) -basic_umap_scatterplot(my_umap$df, #the data frame containing V1 and V2 from UMAP - plot_samples = "some_sample_ID", - colour_by = "DLBCLone_ko") +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/make_alluvial.Rd b/man/make_alluvial.Rd index a956a87..4b41839 100644 --- a/man/make_alluvial.Rd +++ b/man/make_alluvial.Rd @@ -97,7 +97,10 @@ unclassified rates, and flexible labeling and coloring options. } } \examples{ -# Example usage: -# make_alluvial(optimized_result) +\dontrun{ +library(GAMBLR.predict) + +make_alluvial(optimize_params) +} } diff --git a/man/make_neighborhood_plot.Rd b/man/make_neighborhood_plot.Rd index 24b3064..37cabb7 100644 --- a/man/make_neighborhood_plot.Rd +++ b/man/make_neighborhood_plot.Rd @@ -45,6 +45,21 @@ The function extracts the nearest neighbors of the specified sample, draws segme # and 'output' is the result of DLBCLone_predict_single_sample # on sample_id "SAMPLE123": \dontrun{ - make_neighborhood_plot(output, optimization_result$df, "SAMPLE123") +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 index de85d94..f175cab 100644 --- a/man/make_umap_scatterplot.Rd +++ b/man/make_umap_scatterplot.Rd @@ -17,20 +17,39 @@ make_umap_scatterplot( ) } \arguments{ -\item{df}{} +\item{df}{Data frame containing the UMAP coordinates and annotations.} -\item{drop_composite}{} +\item{drop_composite}{If TRUE: removes composite labels from the lymphgen column.} -\item{colour_by}{} +\item{colour_by}{Column name to color points by. Default: "lymphgen".} -\item{drop_other}{} +\item{drop_other}{If TRUE: removes "Other" and "NOS" labels from the lymphgen column.} -\item{high_confidence}{} +\item{high_confidence}{If TRUE: filters the data to include only samples with confidence > 0.7.} -\item{custom_colours}{} +\item{custom_colours}{Custom color palette for the plot. If not provided, uses default GAMBL colors.} -\item{add_labels}{} +\item{add_labels}{If TRUE: adds labels to the points based on the median coordinates of each group.} +} +\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 index 8236c65..bda58ed 100644 --- a/man/nearest_neighbor_heatmap.Rd +++ b/man/nearest_neighbor_heatmap.Rd @@ -40,6 +40,15 @@ based on a DLBCLone model object. Supports outputs from: } } \examples{ -# Assuming 'predicted_out' is a the output of DLBCLone_KNN_predict +# 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/report_accuracy.Rd b/man/report_accuracy.Rd index b6481da..2b5a08a 100644 --- a/man/report_accuracy.Rd +++ b/man/report_accuracy.Rd @@ -45,8 +45,12 @@ Optionally excludes samples assigned to the "Other" class from accuracy calculat } } \examples{ +\dontrun{ +library(GAMBLR.predict) + result <- report_accuracy(predictions_df) result$overall result$per_class +} } From 0e0e1acac2f025186c1c13862dbbbe176b282773 Mon Sep 17 00:00:00 2001 From: Atlantis-Real Date: Wed, 17 Sep 2025 00:02:17 -0700 Subject: [PATCH 3/5] umap.R documentation updates --- NAMESPACE | 1 + R/umap.R | 300 ++++++++++++++++++-------- man/DLBCLone_KNN.Rd | 15 ++ man/DLBCLone_KNN_predict.Rd | 24 ++- man/DLBCLone_ensemble_postprocess.Rd | 14 ++ man/DLBCLone_optimize_params.Rd | 43 ++-- man/DLBCLone_predict.Rd | 35 ++- man/DLBCLone_predict_mixture_model.Rd | 20 +- man/DLBCLone_train_mixture_model.Rd | 16 +- man/assemble_genetic_features.Rd | 30 +++ man/make_and_annotate_umap.Rd | 6 +- man/optimize_purity.Rd | 8 +- man/process_votes.Rd | 7 +- 13 files changed, 395 insertions(+), 124 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 549bb97..d8da5d4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -56,6 +56,7 @@ import(rlang) import(shiny) import(shinybusy) import(shinyjs) +import(stringr) import(tibble) import(tidyr) import(tidyselect) diff --git a/R/umap.R b/R/umap.R index 74ba15c..66e854c 100644 --- a/R/umap.R +++ b/R/umap.R @@ -137,7 +137,6 @@ summarize_all_ssm_status <- function(maf_df, return(mutation_wide) } - #' Assemble genetic features for UMAP input #' #' This function assembles a matrix of genetic features for each sample, including mutation status, @@ -155,7 +154,39 @@ summarize_all_ssm_status <- function(maf_df, #' @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) +#' +#' 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, @@ -316,10 +347,8 @@ assemble_genetic_features <- function(these_samples_metadata, status_combined[,"BCL6_SV"] = 0 status_combined[bcl6_id,"BCL6_SV"] = sv_value return(status_combined) - } - #' Optimize the threshold for classifying samples as "Other" #' #' Performs a post-hoc evaluation of the classification of a sample as one of @@ -345,6 +374,8 @@ assemble_genetic_features <- function(these_samples_metadata, #' (passed to DLBCLone_optimize_params). Default: FALSE #' #' @returns a list of data frames with the predictions and the UMAP input +#' +#' @import dplyr tidyr caret #' @export #' optimize_outgroup <- function(predicted_labels, @@ -413,7 +444,6 @@ optimize_outgroup <- function(predicted_labels, return(best) } - #' Plot the result of a DLBCLone classification #' #' @param test_df Data frame containing the test data with UMAP coordinates @@ -528,7 +558,6 @@ DLBCLone_train_test_plot = function(test_df, 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 @@ -554,12 +583,14 @@ DLBCLone_train_test_plot = function(test_df, #' #' @examples #' -#' #library(GAMBLR.predict) +#' 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")) -#' my_umap <- make_and_annotate_umap( +#' +#' make_umap <- make_and_annotate_umap( #' df=all_full_status, #' metadata=dlbcl_meta #' ) @@ -614,7 +645,6 @@ make_and_annotate_umap = function(df, } } - if(missing(df)){ stop("provide a data frame or matrix with one row for each sample and a numeric column for each mutation feature") @@ -748,8 +778,7 @@ make_and_annotate_umap = function(df, rng_type = "deterministic" ) #IMPORTANT: n_threads must never be changed because it will break reproducibility - } - + } }else{ message("transforming each data point individually using the provided UMAP model. This will take some time.") @@ -771,7 +800,6 @@ make_and_annotate_umap = function(df, } - if(model_provided){ # model was generated here message("model given to function") @@ -785,8 +813,6 @@ make_and_annotate_umap = function(df, umap_df = left_join(umap_df,metadata,by=join_column) } results = list() - - results[["df"]]=umap_df results[["features"]] = df @@ -835,11 +861,16 @@ make_and_annotate_umap = function(df, #' - 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. #' -#' @examples -#' # Example usage: -#' # result <- process_votes(knn_output_df, k = 7) -#' +#' @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"), @@ -971,8 +1002,6 @@ process_votes <- function(df, } - - #' 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. @@ -990,12 +1019,17 @@ process_votes <- function(df, #' - 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 caret -#' @examples -#' # Example usage: -#' # result <- optimize_purity(processed_votes, prediction_column = "pred_label", truth_column = "true_label") -#' +#' @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, @@ -1147,7 +1181,6 @@ optimize_purity <- function(optimized_model_object, 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. @@ -1172,17 +1205,34 @@ optimize_purity <- function(optimized_model_object, #' - 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 = train_features, -#' test_df = new_samples, -#' metadata = sample_metadata, -#' DLBCLone_KNN_out = model_out +#' train_df = all_full_status, +#' test_df = lyseq_validation_data, +#' metadata = dlbcl_meta, +#' DLBCLone_KNN_out = dlbcl_knn, +#' mode = "batch" +#' ) #' -#' @export #' DLBCLone_KNN_predict <- function(train_df, test_df, @@ -1243,7 +1293,6 @@ DLBCLone_KNN_predict <- function(train_df, } } - #' Run DLBCLone KNN Classification #' #' Weighted KNN on a feature (mutation) matrix with optional upweighting of @@ -1290,7 +1339,23 @@ DLBCLone_KNN_predict <- function(train_df, #' \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, @@ -1886,8 +1951,6 @@ DLBCLone_KNN <- function(features_df, 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 @@ -1925,34 +1988,47 @@ DLBCLone_KNN <- function(features_df, #' 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 caret -#' +#' @import dplyr tidyr 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 +#' ) #' #' # Aim to maximize classification of samples into non-Other class -#' lymphgen_lyseq_no_other = -#' GAMBLR.predict::DLBCLone_optimize_params( -#' dlbcl_status_combined_lyseq, -#' dlbcl_meta_lyseq_train, -#' min_k = 5,max_k=23, -#' optimize_for_other = F, -#' truth_classes = c("MCD","EZB","ST2","BN2")) +#' 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") #' -#' DLBCLone_lymphgen_lyseq_prime_opt = -#' GAMBLR.predict::DLBCLone_optimize_params( -#' dlbcl_status_combined_lyseq_prime, -#' dlbcl_meta_lyseq_train, -#' min_k = 5,max_k=23, -#' optimize_for_other = T, -#' truth_classes = c("MCD","EZB","ST2","BN2","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, @@ -2055,7 +2131,6 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, } } - for(use_w in weights_opt){ for(k in ks){ best_w_acc = NULL @@ -2068,8 +2143,6 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, rownames(train_coords) = train_ids rownames(test_coords) = test_ids - - pred = weighted_knn_predict_with_conf( train_coords = train_coords, train_labels = train_labels, @@ -2189,7 +2262,7 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, if(verbose){ print("true_factor") - print(levels(true_factor )) + print(levels(true_factor)) print("===") print(levels(xx_d$predicted_label)) } @@ -2324,7 +2397,6 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, 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, @@ -2348,7 +2420,6 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, by_score, "Other")) - if(optimize_for_other){ pred = mutate(pred,predicted_label_optimized = ifelse(other_score > best_params$threshold_outgroup, @@ -2433,14 +2504,23 @@ DLBCLone_optimize_params = function(combined_mutation_status_df, #' - predictions: updated predictions_df (DLBCLone_wo or DLBCLone_wc / DLBCLone_wc_simplified) #' - consistency_report: per-sample counts and decision flags #' -#' @import dplyr tidyr +#' @import dplyr tidyr tibble #' @export +#' +#' @examples +#' \dontrun{ +#' library(GAMBLR.predict) +#' +#' DLBCLone_ensemble_postprocess( +#' optimized_model = optimize_params_T, +#' other_min = 3, +#' assign_composites = TRUE, +#' any_split = TRUE, +#' optimize_per_class = TRUE +#' ) +#' } +#' DLBCLone_ensemble_postprocess <- function( - #all_predictions, - #truth_classes, - #predictions_df, - #outs_df, - #truth_column, optimized_model, assign_composites = FALSE, other_min = 2, @@ -2626,9 +2706,6 @@ DLBCLone_ensemble_postprocess <- function( return(to_return) } - - - #' Weighted k-nearest neighbor with confidence estimate #' #' @param train_coords Data frame of coordinates for labeled (training) samples. @@ -2672,7 +2749,7 @@ DLBCLone_ensemble_postprocess <- function( #' \item{total_w}{sum of weights for in-group neighbors used} #' \item{pred_w}{weight supporting the predicted class} #' -#' @import FNN +#' @import FNN dplyr tibble tidyr #' @export #' #' @@ -2937,18 +3014,37 @@ prepare_single_sample_DLBCLone <- function(optimized_model,seed=12345){ #' samples to estimate overall accuracy #' @param truth_classes Vector of classes to use for training and testing. Default: c("EZB","MCD","ST2","N1","BN2") #' @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 -#' predict_single_sample_DLBCLone( -#' seed = 1234, -#' test_df = test_df, -#' train_df = train_df, -#' train_metadata = train_metadata, -#' umap_out = umap_out, -#' best_params = best_params -#' predictions_df = predictions_df, -#' annotate_accuracy = TRUE +#' +#' 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 <- DLBCLone_predict( +#' mutation_status = optimize_params$features[1,], +#' optimized_model = optimize_params #' ) #' DLBCLone_predict <- function( @@ -3097,8 +3193,6 @@ DLBCLone_predict <- function( }else{ stop("mode must be one of DLBCLone_w or DLBCLone_i") } - - # coords for outputs (always from the unified projection) predictions_test_df <- dplyr::left_join( @@ -3164,7 +3258,6 @@ DLBCLone_predict <- function( return(to_return) } - #' 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". @@ -3186,12 +3279,26 @@ DLBCLone_predict <- function( #' \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 -#' result <- DLBCLone_train_mixture_model(umap_out) +#' +#' 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) -#' @import mclust -#' -#' @export +#' DLBCLone_train_mixture_model = function(umap_out, probability_threshold = 0.5, density_max_threshold = 0.05, @@ -3293,13 +3400,32 @@ return(to_return) #' \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: -#' result <- DLBCLone_predict_mixture_model(model, umap_out) +#' +#' 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$gaussian_mixture_model, +#' umap_out=make_umap +#' ) #' head(result$predictions) #' -#' @import mclust -#' @export DLBCLone_predict_mixture_model = function(model, umap_out, probability_threshold = 0.5, @@ -3338,17 +3464,14 @@ 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, @@ -3356,4 +3479,3 @@ to_return = list(gaussian_mixture_model = gmm_supervised, ) return(to_return) } - diff --git a/man/DLBCLone_KNN.Rd b/man/DLBCLone_KNN.Rd index 531a8d6..76c2048 100644 --- a/man/DLBCLone_KNN.Rd +++ b/man/DLBCLone_KNN.Rd @@ -89,3 +89,18 @@ 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 index ccfc263..1f7b9a6 100644 --- a/man/DLBCLone_KNN_predict.Rd +++ b/man/DLBCLone_KNN_predict.Rd @@ -52,10 +52,26 @@ for more stable results when predicting multiple samples. # 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 = train_features, - test_df = new_samples, - metadata = sample_metadata, - DLBCLone_KNN_out = model_out + 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_ensemble_postprocess.Rd b/man/DLBCLone_ensemble_postprocess.Rd index caa720d..e52deda 100644 --- a/man/DLBCLone_ensemble_postprocess.Rd +++ b/man/DLBCLone_ensemble_postprocess.Rd @@ -44,3 +44,17 @@ Post-process KNN results across K to score consistency, (optionally) refine classified/Other cutoffs per-class and (optionally) assign composite classes } +\examples{ +\dontrun{ +library(GAMBLR.predict) + +DLBCLone_ensemble_postprocess( + optimized_model = optimize_params_T, + other_min = 3, + assign_composites = TRUE, + any_split = TRUE, + optimize_per_class = TRUE +) +} + +} diff --git a/man/DLBCLone_optimize_params.Rd b/man/DLBCLone_optimize_params.Rd index 02cffed..ea91618 100644 --- a/man/DLBCLone_optimize_params.Rd +++ b/man/DLBCLone_optimize_params.Rd @@ -73,25 +73,40 @@ Optimize parameters for classifying samples using UMAP and k-nearest neighbor } \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 +) # Aim to maximize classification of samples into non-Other class -lymphgen_lyseq_no_other = -GAMBLR.predict::DLBCLone_optimize_params( - dlbcl_status_combined_lyseq, - dlbcl_meta_lyseq_train, - min_k = 5,max_k=23, - optimize_for_other = F, - truth_classes = c("MCD","EZB","ST2","BN2")) +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") -DLBCLone_lymphgen_lyseq_prime_opt = -GAMBLR.predict::DLBCLone_optimize_params( - dlbcl_status_combined_lyseq_prime, - dlbcl_meta_lyseq_train, - min_k = 5,max_k=23, - optimize_for_other = T, - truth_classes = c("MCD","EZB","ST2","BN2","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.Rd b/man/DLBCLone_predict.Rd index f69236e..c5f870d 100644 --- a/man/DLBCLone_predict.Rd +++ b/man/DLBCLone_predict.Rd @@ -45,15 +45,32 @@ a list of data frames with the predictions, the UMAP input, the model, and a ggp Predict class for one or more samples using a pre-trained DLBCLone model } \examples{ -predict_single_sample_DLBCLone( - seed = 1234, - test_df = test_df, - train_df = train_df, - train_metadata = train_metadata, - umap_out = umap_out, - best_params = best_params - predictions_df = predictions_df, - annotate_accuracy = TRUE + +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 <- DLBCLone_predict( + mutation_status = optimize_params$features[1,], + optimized_model = optimize_params ) } diff --git a/man/DLBCLone_predict_mixture_model.Rd b/man/DLBCLone_predict_mixture_model.Rd index 3d94ea5..198bb61 100644 --- a/man/DLBCLone_predict_mixture_model.Rd +++ b/man/DLBCLone_predict_mixture_model.Rd @@ -45,7 +45,25 @@ Assigns class predictions and optionally reclassifies samples as "Other" based o } \examples{ # Predict on new UMAP data using a trained mixture model: -result <- DLBCLone_predict_mixture_model(model, umap_out) + +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$gaussian_mixture_model, + umap_out=make_umap +) head(result$predictions) } diff --git a/man/DLBCLone_train_mixture_model.Rd b/man/DLBCLone_train_mixture_model.Rd index 4a83119..454382e 100644 --- a/man/DLBCLone_train_mixture_model.Rd +++ b/man/DLBCLone_train_mixture_model.Rd @@ -41,6 +41,20 @@ Assigns class predictions and optionally reclassifies samples as "Other" based o } } \examples{ -result <- DLBCLone_train_mixture_model(umap_out) + +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/assemble_genetic_features.Rd b/man/assemble_genetic_features.Rd index 4dbc36b..292b57d 100644 --- a/man/assemble_genetic_features.Rd +++ b/man/assemble_genetic_features.Rd @@ -50,3 +50,33 @@ Matrix of assembled features for each sample. 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/make_and_annotate_umap.Rd b/man/make_and_annotate_umap.Rd index 3fbdf52..777098e 100644 --- a/man/make_and_annotate_umap.Rd +++ b/man/make_and_annotate_umap.Rd @@ -64,12 +64,14 @@ Run UMAP and attach result to metadata } \examples{ -#library(GAMBLR.predict) +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")) -my_umap <- make_and_annotate_umap( + +make_umap <- make_and_annotate_umap( df=all_full_status, metadata=dlbcl_meta ) diff --git a/man/optimize_purity.Rd b/man/optimize_purity.Rd index 03170ad..221550b 100644 --- a/man/optimize_purity.Rd +++ b/man/optimize_purity.Rd @@ -41,7 +41,11 @@ The function returns the best accuracy achieved and the corresponding purity thr } } \examples{ -# Example usage: -# result <- optimize_purity(processed_votes, prediction_column = "pred_label", truth_column = "true_label") + +\dontrun{ +library(GAMBLR.predict) + +result <- optimize_purity(processed_votes, prediction_column = "pred_label", truth_column = "true_label") +} } diff --git a/man/process_votes.Rd b/man/process_votes.Rd index 050ec5d..6583f1a 100644 --- a/man/process_votes.Rd +++ b/man/process_votes.Rd @@ -52,7 +52,10 @@ The function also supports custom logic for handling the "Other" class, includin } } \examples{ -# Example usage: -# result <- process_votes(knn_output_df, k = 7) +/dontrun{ +library(GAMBLR.predict) + +result <- process_votes(knn_output_df, k = 7) +} } From 2baa8698a9c0777c5adf792dd36df6dbd73adb2b Mon Sep 17 00:00:00 2001 From: Atlantis-Real Date: Wed, 17 Sep 2025 20:37:00 -0700 Subject: [PATCH 4/5] removing summarize_all_ssm_status, updating assemble_genetic_features and fixing its SV spacing issue --- NAMESPACE | 1 - R/umap.R | 250 +++++++------------------------ man/assemble_genetic_features.Rd | 43 ++---- man/summarize_all_ssm_status.Rd | 83 ---------- 4 files changed, 73 insertions(+), 304 deletions(-) delete mode 100644 man/summarize_all_ssm_status.Rd diff --git a/NAMESPACE b/NAMESPACE index d8da5d4..2fadee4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -31,7 +31,6 @@ 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) diff --git a/R/umap.R b/R/umap.R index 66e854c..3e1f37c 100644 --- a/R/umap.R +++ b/R/umap.R @@ -1,142 +1,3 @@ -#' Summarize SSM (Somatic Single Nucleotide Mutation) Status Across Samples -#' -#' 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. -#' -#' @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 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 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 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. -#' -#' @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{ -#' 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 -#' ) -#' } -#' -#' @import dplyr tidyr tibble readr -#' @export -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, @@ -159,36 +20,28 @@ summarize_all_ssm_status <- function(maf_df, #' @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) -#' -#' 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 -#' ) -#' } -#' +#' all_meta = get_gambl_metadata() %>% +#' dplyr::filter(pathology=="DLBCL",seq_type=="genome") +#' +#' all_maf = get_all_coding_ssm(all_meta,include_silent=TRUE) +#' +#' sv_all =get_combined_sv(all_meta) +#' +#' anno_sv = annotate_sv(sv_all) +#' +#' feat_mat = assemble_genetic_features(all_meta, +#' genes=c("EZH2","SOCS1","PIM1", +#' "MYD88","CREBBP","SGK1", +#' "NOTCH2","NOTCH1"), +#' maf_with_synon = all_maf, +#' synon_genes=c("PIM1","SOCS1")) +#'} +#' @export assemble_genetic_features <- function(these_samples_metadata, - metadata_columns = c("bcl2_ba","bcl6_ba","myc_ba"), + sv_from_metadata = c(BCL2="bcl2_ba", + BCL6="bcl6_ba", + MYC="myc_ba"), genes, synon_genes, maf_with_synon, @@ -197,7 +50,7 @@ assemble_genetic_features <- function(these_samples_metadata, sv_value = 3, synon_value = 1, coding_value = 2, - include_ashm = TRUE, + include_ashm = FALSE, annotated_sv, include_GAMBL_sv= TRUE, review_hotspots = TRUE, @@ -317,38 +170,49 @@ assemble_genetic_features <- function(these_samples_metadata, 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) - + sv_samples = list() + if(!is.null(sv_from_metadata)){ + for(oncogene in names(sv_from_metadata)){ + metadata_column = sv_from_metadata[[oncogene]] + #assume POS/NEG encoding. Warn if not consistent with this + col_vals = unique(these_samples_metadata[[metadata_column]]) + if(!"POS" %in% col_vals){ + stop(paste("column",metadata_column, "doesn't have any POS values", + "SV encoding is expected to be POS/NEG")) + } + sv_samples[[oncogene]] = these_samples_metadata[which(these_samples_metadata[[metadata_column]]=="POS"),] %>% pull(sample_id) + + } + } + #print(sv_samples) + for(oncogene in names(sv_from_metadata)){ + message(paste("SVs in",oncogene,"from metadata:",length(sv_samples[[oncogene]]))) + } + # 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)) + for(oncogene in names(sv_from_metadata)){ + print(paste("oncogene:",oncogene)) + sv_id = filter(annotated_sv,!is.na(partner),gene==oncogene) %>% pull(tumour_sample_id) + if(oncogene %in% names(sv_samples)){ + sv_samples[[oncogene]] = union( sv_samples[[oncogene]],sv_id) + } + } + for(oncogene in names(sv_from_metadata)){ + message(paste("SVs in",oncogene,"after adding annotated SVs:",length(sv_samples[[oncogene]]))) + } + } + for(oncogene in names(sv_samples)){ + onco_column = paste(oncogene,"_SV", sep = "") + status_combined[[onco_column]] = 0 + status_combined[sv_samples[[oncogene]],onco_column] = sv_value } - - - 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) + } + #' Optimize the threshold for classifying samples as "Other" #' #' Performs a post-hoc evaluation of the classification of a sample as one of diff --git a/man/assemble_genetic_features.Rd b/man/assemble_genetic_features.Rd index 292b57d..0461c28 100644 --- a/man/assemble_genetic_features.Rd +++ b/man/assemble_genetic_features.Rd @@ -6,7 +6,7 @@ \usage{ assemble_genetic_features( these_samples_metadata, - metadata_columns = c("bcl2_ba", "bcl6_ba", "myc_ba"), + sv_from_metadata = c(BCL2 = "bcl2_ba", BCL6 = "bcl6_ba", MYC = "myc_ba"), genes, synon_genes, maf_with_synon, @@ -15,7 +15,7 @@ assemble_genetic_features( sv_value = 3, synon_value = 1, coding_value = 2, - include_ashm = TRUE, + include_ashm = FALSE, annotated_sv, include_GAMBL_sv = TRUE, review_hotspots = TRUE, @@ -25,8 +25,6 @@ assemble_genetic_features( \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.} @@ -42,6 +40,8 @@ assemble_genetic_features( \item{coding_value}{Value to assign for coding mutations (default: 2).} \item{verbose}{Defaults to FALSE} + +\item{metadata_columns}{Columns in metadata to use for SV status (default: c("bcl2_ba","bcl6_ba","myc_ba")).} } \value{ Matrix of assembled features for each sample. @@ -51,32 +51,21 @@ This function assembles a matrix of genetic features for each sample, including aSHM counts, and structural variant status for BCL2, BCL6, and MYC. It supports both genome and capture sequencing types. } \examples{ - \dontrun{ -library(GAMBLR.predict) +all_meta = get_gambl_metadata() \%>\% + dplyr::filter(pathology=="DLBCL",seq_type=="genome") -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_maf = get_all_coding_ssm(all_meta,include_silent=TRUE) -all_sv_anno = get_combined_sv(these_samples_metadata = dlbcl_meta) +sv_all =get_combined_sv(all_meta) -all_sv_anno = annotate_sv(all_sv_anno) +anno_sv = annotate_sv(sv_all) -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 -) -} - +feat_mat = assemble_genetic_features(all_meta, + genes=c("EZH2","SOCS1","PIM1", + "MYD88","CREBBP","SGK1", + "NOTCH2","NOTCH1"), + maf_with_synon = all_maf, + synon_genes=c("PIM1","SOCS1")) +} } diff --git a/man/summarize_all_ssm_status.Rd b/man/summarize_all_ssm_status.Rd deleted file mode 100644 index ec7b844..0000000 --- a/man/summarize_all_ssm_status.Rd +++ /dev/null @@ -1,83 +0,0 @@ -% 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{ -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 -) -} - -} From 844cf41d9566638228cbe9384a22200344f3eaee Mon Sep 17 00:00:00 2001 From: Atlantis-Real Date: Thu, 18 Sep 2025 15:39:55 -0700 Subject: [PATCH 5/5] updates to DLBCLone_predict and nearest_neighbor_heatmap with hashed comments --- R/plotting.R | 9 ++++++--- R/umap.R | 1 + 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/R/plotting.R b/R/plotting.R index 4ae2b84..728d073 100644 --- a/R/plotting.R +++ b/R/plotting.R @@ -168,11 +168,14 @@ nearest_neighbor_heatmap <- function( if ("metadata" %in% names(DLBCLone_model) && !is.null(DLBCLone_model$metadata)) { if(DLBCLone_model$type=="DLBCLone_predict"){ #need to pool together training sample metadata with predictions - train_meta = DLBCLone_model$metadata %>% - dplyr::select(dplyr::all_of(c("sample_id", truth_column, metadata_cols))) + train_meta = DLBCLone_model$optimized_predictions %>% ############################################ metadata -> optimized_predictions + dplyr::select(dplyr::all_of(c("sample_id", truth_column, pred_col, metadata_cols))) #################################################### pred_col added #preds <- dplyr::left_join(preds, DLBCLone_model$metadata, by = "sample_id") - preds <- dplyr::bind_rows( + colnames(preds) + colnames(train_meta) + + preds <- dplyr::left_join( ######################################################## bind_rows -> left_join train_meta, preds ) %>% dplyr::distinct(sample_id, .keep_all = TRUE) diff --git a/R/umap.R b/R/umap.R index 3e1f37c..b337a31 100644 --- a/R/umap.R +++ b/R/umap.R @@ -3112,6 +3112,7 @@ DLBCLone_predict <- function( #anno_df = dplyr::left_join(predictions_df, train_metadata, by = "sample_id"), anno_out = anno_out, type = "DLBCLone_predict", + optimized_predictions = optimized_model$predictions, # <- NECESSARY FOR HEATMAP!! unprocessed_votes = unprocessed ) if(weight_core_features){