From 01e7f5b83a93d6a92c5c3c04d8fa48f351fc5e5e Mon Sep 17 00:00:00 2001 From: Konstantinos Antonopoulos Date: Tue, 9 Dec 2025 14:39:35 +0100 Subject: [PATCH 1/4] Update version from 1.0.0 to 1.0.1 --- DESCRIPTION | 2 +- README.md | 2 +- index.md | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d1e9fba..ce748c4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: HDAnalyzeR Title: Streamlining Data Analysis for Biomarker Research -Version: 1.0.0 +Version: 1.0.1 Authors@R: c( person("Konstantinos", "Antonopoulos", , "k.antono@outlook.com", diff --git a/README.md b/README.md index e10748d..b7db870 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [![R-CMD-check](https://github.com/kantonopoulos/HDAnalyzeR/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/kantonopoulos/HDAnalyzeR/actions/workflows/R-CMD-check.yaml) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.13786124.svg)](https://doi.org/10.5281/zenodo.13786124) -[![Version](https://img.shields.io/badge/Version-1.0.0-purple)](https://github.com/kantonopoulos/HDAnalyzeR) +[![Version](https://img.shields.io/badge/Version-1.0.1-purple)](https://github.com/kantonopoulos/HDAnalyzeR) [![License](https://img.shields.io/badge/license-Apache2.0-yellow)](https://github.com/kantonopoulos/HDAnalyzeR/blob/main/LICENSE.md) HDAnalyzeR is an R package developed by the Human Disease Blood Atlas to simplify data analysis for biomarker discovery in disease research. It provides a set of user-friendly functions to efficiently process omics data and gain insights into disease mechanisms. With HDAnalyzeR, you can easily perform tasks like protein differential expression analysis, classification modeling, dimensionality reduction, imputation, and data visualization. diff --git a/index.md b/index.md index 9520d72..cc38499 100644 --- a/index.md +++ b/index.md @@ -2,7 +2,7 @@ [![R-CMD-check](https://github.com/kantonopoulos/HDAnalyzeR/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/kantonopoulos/HDAnalyzeR/actions/workflows/R-CMD-check.yaml) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.13786124.svg)](https://doi.org/10.5281/zenodo.13786124) -[![Version](https://img.shields.io/badge/Version-1.0.0-purple)](https://github.com/kantonopoulos/HDAnalyzeR) +[![Version](https://img.shields.io/badge/Version-1.0.1-purple)](https://github.com/kantonopoulos/HDAnalyzeR) [![License](https://img.shields.io/badge/license-Apache2.0-yellow)](https://github.com/kantonopoulos/HDAnalyzeR/blob/main/LICENSE.md) HDAnalyzeR is an R package developed by the Human Disease Blood Atlas to simplify data analysis for biomarker discovery in disease research. It provides a set of user-friendly functions to efficiently process omics data and gain insights into disease mechanisms. With HDAnalyzeR, you can easily perform tasks like protein differential expression analysis, classification modeling, dimensionality reduction, imputation, and data visualization. From 842ec28039aa417e3f33873435659b33fad19a4b Mon Sep 17 00:00:00 2001 From: Konstantinos Antonopoulos Date: Tue, 9 Dec 2025 16:22:59 +0100 Subject: [PATCH 2/4] Add feature importance for each class for glmnet multiclass models --- .Rbuildignore | 2 + R/classification_models.R | 136 +++- R/visualize.R | 8 +- inst/doc/HDAnalyzeR.html | 813 ---------------------- inst/doc/before_using_hdanalyzer.html | 461 ------------- inst/doc/classification.html | 926 -------------------------- inst/doc/data_imputation.html | 628 ----------------- inst/doc/data_preprocessing.html | 590 ---------------- inst/doc/differential_expression.html | 568 ---------------- inst/doc/post_analysis.html | 540 --------------- man/hd_plot_feature_heatmap.Rd | 2 +- man/variable_imp.Rd | 3 + 12 files changed, 118 insertions(+), 4559 deletions(-) delete mode 100644 inst/doc/HDAnalyzeR.html delete mode 100644 inst/doc/before_using_hdanalyzer.html delete mode 100644 inst/doc/classification.html delete mode 100644 inst/doc/data_imputation.html delete mode 100644 inst/doc/data_preprocessing.html delete mode 100644 inst/doc/differential_expression.html delete mode 100644 inst/doc/post_analysis.html diff --git a/.Rbuildignore b/.Rbuildignore index 3ab49c4..008e79b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -17,3 +17,5 @@ ^index.md ^case_studies$ ^vignettes$ +^docs$ +^pkgdown$ diff --git a/R/classification_models.R b/R/classification_models.R index 5b5bef5..18073f2 100644 --- a/R/classification_models.R +++ b/R/classification_models.R @@ -1003,13 +1003,18 @@ generate_title <- function(features, if ("features" %in% title) { title_parts <- c(title_parts, paste0('Features = ', - nrow(features |> dplyr::filter(!!rlang::sym("Scaled_Importance") > 0)), + nrow(features |> + dplyr::filter(!!rlang::sym("Scaled_Importance") > 0) |> + dplyr::select(dplyr::any_of(c("Feature"))) |> + dplyr::distinct()), ' ')) } if ("top-features" %in% title) { title_parts <- c(title_parts, paste0('Top-features = ', - nrow(features |> dplyr::filter(!!rlang::sym("Scaled_Importance") >= 0.5)), + nrow(features |> dplyr::filter(!!rlang::sym("Scaled_Importance") >= 0.5) |> + dplyr::select(dplyr::any_of(c("Feature"))) |> + dplyr::distinct()), ' ')) } @@ -1034,6 +1039,7 @@ generate_title <- function(features, #' @param mixture The mixture parameter for the elastic net. If NULL it will be tuned. Default is NULL. #' @param palette The color palette for the classes. If it is a character, it should be one of the palettes from `hd_palettes()`. Default is NULL. #' @param title Vector of title elements to include in the plot. +#' @param engine Either glmnet, rf or lr. If glmnet and the model type is multi-class then the variable importance is calculated per class. #' @param seed Seed for reproducibility. Default is 123. #' #' @return A model object containing the features and the feature importance plot. @@ -1051,6 +1057,7 @@ variable_imp <- function(dat, "features", "top-features"), verbose = TRUE, + engine = 'glmnet', seed = 123) { if (verbose){ @@ -1062,19 +1069,48 @@ variable_imp <- function(dat, metrics <- dat[["metrics"]] mixture <- dat[["mixture"]] model_type <- dat[["model_type"]] - - features <- final |> - workflows::extract_fit_parsnip() |> - vip::vi() |> - # Clip negatives to zero to avoid negative importance in Random Forest permutation testing - dplyr::mutate(Importance = dplyr::if_else(!!rlang::sym("Importance") < 0, - 0, - !!rlang::sym("Importance")), - Variable = forcats::fct_reorder(Variable, !!rlang::sym("Importance"))) |> - dplyr::arrange(dplyr::desc(!!rlang::sym("Importance"))) |> - # Min max scaling with min = 0 always and max = 1 - dplyr::mutate(Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance"))) |> - dplyr::rename(Feature = !!rlang::sym("Variable")) + + if (engine == 'glmnet') { + features <- final |> + workflows::extract_fit_parsnip() |> + broom::tidy() |> + dplyr::filter(!!rlang::sym("term") != "(Intercept)") |> + dplyr::select(-dplyr::any_of(c("penalty"))) |> + dplyr::mutate( + Feature = !!rlang::sym("term"), + Importance = abs(!!rlang::sym("estimate")) + ) |> + dplyr::arrange(dplyr::desc(!!rlang::sym("Importance"))) + + # Scale importance + if (model_type == "multi_class") { + features <- features |> + dplyr::group_by(class) |> + dplyr::mutate( + Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance")), + Feature_plot = tidytext::reorder_within(!!rlang::sym("Feature"), !!rlang::sym("Importance"), !!rlang::sym("class")) + ) |> + dplyr::ungroup() + } else { + features <- features |> + dplyr::mutate(Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance"))) |> + dplyr::mutate(Feature = forcats::fct_reorder(!!rlang::sym("Feature"), !!rlang::sym("Importance"))) + } + + } else { + features <- final |> + workflows::extract_fit_parsnip() |> + vip::vi() |> + # Clip negatives to zero to avoid negative importance in Random Forest permutation testing + dplyr::mutate(Importance = dplyr::if_else(!!rlang::sym("Importance") < 0, + 0, + !!rlang::sym("Importance")), + Variable = forcats::fct_reorder(!!rlang::sym("Variable"), !!rlang::sym("Importance"))) |> + dplyr::arrange(dplyr::desc(!!rlang::sym("Importance"))) |> + # Min max scaling with min = 0 always and max = 1 + dplyr::mutate(Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance"))) |> + dplyr::rename(Feature = !!rlang::sym("Variable")) + } if (model_type == "binary_class") { @@ -1096,6 +1132,27 @@ variable_imp <- function(dat, pal <- c("#883268") } + } else if (model_type == "multi_class" && engine == "glmnet") { + + title_text <- generate_title(features = features, + accuracy = as.numeric(metrics[["accuracy"]]), + sensitivity = as.numeric(metrics[["sensitivity"]]), + specificity = as.numeric(metrics[["specificity"]]), + auc = NULL, + mixture = mixture, + title = title) + + pals <- hd_palettes() + if (!is.null(palette) && is.null(names(palette))) { + pal <- pals[palette] + pal <- unlist(pals[[palette]]) + } else if (!is.null(palette)) { + pal <- palette + } else { + classes <- unique(features$class) + pal <- rep("#883268", length(classes)) + } + } else if (model_type == "multi_class") { title_text <- generate_title(features = features, @@ -1118,18 +1175,36 @@ variable_imp <- function(dat, pal <- c("#883268") case <- "case" } - - var_imp_plot <- features |> - dplyr::filter(!!rlang::sym("Scaled_Importance") > 0) |> - ggplot2::ggplot(ggplot2::aes(x = !!rlang::sym("Scaled_Importance"), y = !!rlang::sym("Feature"))) + - ggplot2::geom_col(ggplot2::aes(fill = ifelse(!!rlang::sym("Scaled_Importance") > 0.5, case, NA))) + - ggplot2::labs(y = NULL) + - ggplot2::scale_x_continuous(breaks = c(0, 1), expand = c(0, 0)) + # Keep x-axis tick labels at 0 and 1 - ggplot2::scale_fill_manual(values = pal, na.value = "grey80") + - ggplot2::ggtitle(label = title_text) + - ggplot2::xlab('Importance') + - ggplot2::ylab('Features') + - theme_hd() + + if (model_type == "multi_class" && engine == "glmnet") { + var_imp_plot <- features |> + dplyr::filter(!!rlang::sym("Scaled_Importance") > 0) |> + ggplot2::ggplot(ggplot2::aes(x = !!rlang::sym("Feature_plot"), y = !!rlang::sym("Scaled_Importance"))) + + ggplot2::geom_col(ggplot2::aes(fill = ifelse(!!rlang::sym("Scaled_Importance") > 0.5, class, NA))) + + ggplot2::facet_wrap(~ class, scales = "free_y") + + tidytext::scale_x_reordered() + # cleans axis labels + ggplot2::coord_flip() + # horizontal bars + ggplot2::scale_fill_manual(values = pal, na.value = "grey80") + + ggplot2::scale_y_continuous(expand = c(0,0), limits = c(0,1)) + + ggplot2::labs(x = "Feature", y = "Importance", title = title_text) + + theme_hd(angled = 90) + + ggplot2::theme( + axis.text.x = ggplot2::element_text(hjust = 0.5), + panel.spacing = ggplot2::unit(1, "lines") + ) + } else { + var_imp_plot <- features |> + dplyr::filter(!!rlang::sym("Scaled_Importance") > 0) |> + ggplot2::ggplot(ggplot2::aes(x = !!rlang::sym("Scaled_Importance"), y = !!rlang::sym("Feature"))) + + ggplot2::geom_col(ggplot2::aes(fill = ifelse(!!rlang::sym("Scaled_Importance") > 0.5, case, NA))) + + ggplot2::labs(y = NULL) + + ggplot2::scale_x_continuous(breaks = c(0, 1), expand = c(0, 0)) + # Keep x-axis tick labels at 0 and 1 + ggplot2::scale_fill_manual(values = pal, na.value = "grey80") + + ggplot2::ggtitle(label = title_text) + + ggplot2::xlab('Importance') + + ggplot2::ylab('Features') + + theme_hd() + } if (isFALSE(y_labels)) { var_imp_plot <- var_imp_plot + @@ -1297,6 +1372,7 @@ hd_model_rreg <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'glmnet', seed = seed) } else if (dat[["model_type"]] == "multi_class") { @@ -1315,6 +1391,7 @@ hd_model_rreg <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'glmnet', seed = seed) } else { @@ -1333,6 +1410,7 @@ hd_model_rreg <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'glmnet', seed = seed) } @@ -1478,6 +1556,7 @@ hd_model_rf <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'rf', seed = seed) } else if (dat[["model_type"]] == "multi_class") { dat <- evaluate_multiclass_model(dat = dat, @@ -1494,6 +1573,7 @@ hd_model_rf <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'rf', seed = seed) } else { dat <- evaluate_regression_model(dat = dat, @@ -1510,6 +1590,7 @@ hd_model_rf <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'rf', seed = seed) } @@ -1641,6 +1722,7 @@ hd_model_lr <- function(dat, y_labels = plot_y_labels, title = plot_title, verbose = verbose, + engine = 'lr', seed = seed) if (dat[["features"]] |> nrow() < 3) { diff --git a/R/visualize.R b/R/visualize.R index c090456..e6446f7 100644 --- a/R/visualize.R +++ b/R/visualize.R @@ -320,7 +320,7 @@ hd_plot_regression <- function(dat, #' `hd_plot_feature_heatmap` plots a summary bubble-heatmap of the combined differential #' expression and classification model results. The heatmap shows the log2 fold change #' and adjusted p-value of the differential expression results, and the scaled importance -#' and sign of the classification model results. The heatmap is ordered and the selected +#' of the classification model results. The heatmap is ordered and the selected #' assays are based on the specified control group. * are added to the features that pass #' the p-value threshold. #' @@ -414,8 +414,7 @@ hd_plot_feature_heatmap <- function(de_results, res_model <- model_results[[i]][["features"]] |> dplyr::filter(!!rlang::sym("Feature") %in% assays) |> dplyr::select(!!rlang::sym("Feature"), - !!rlang::sym("Scaled_Importance"), - !!rlang::sym("Sign")) + !!rlang::sym("Scaled_Importance")) res_combined <- res_de |> dplyr::left_join(res_model, by = c("Feature")) |> @@ -434,11 +433,10 @@ hd_plot_feature_heatmap <- function(de_results, ) |> ggplot2::ggplot(ggplot2::aes(x = !!rlang::sym("Feature"), y = !!rlang::sym("control_group"))) + ggplot2::geom_tile(ggplot2::aes(fill = !!rlang::sym("logFC")), color = "white") + - ggplot2::geom_point(ggplot2::aes(size = !!rlang::sym("Scaled_Importance"), color = !!rlang::sym("Sign"))) + + ggplot2::geom_point(ggplot2::aes(size = !!rlang::sym("Scaled_Importance"))) + ggplot2::geom_point(ggplot2::aes(size = !!rlang::sym("Scaled_Importance")), shape = 1, colour = "black") + ggplot2::geom_text(ggplot2::aes(label = ifelse(!!rlang::sym("adj.P.Val") < pval_lim, "*", "")), color = "black", size = 3) + ggplot2::scale_fill_gradient2(low = "#317EC2", mid = "white", high = "#C03830", midpoint = 0, name = "Log2 FC") + - ggplot2::scale_color_manual(values = c("NEG" = "#317EC2", "POS" = "#C03830"), name = "Sign", na.translate = FALSE) + ggplot2::scale_size(name = "Importance") + ggplot2::labs(x = "Feature", y = "Control Group") + ggplot2::theme_minimal() + diff --git a/inst/doc/HDAnalyzeR.html b/inst/doc/HDAnalyzeR.html deleted file mode 100644 index e2192d6..0000000 --- a/inst/doc/HDAnalyzeR.html +++ /dev/null @@ -1,813 +0,0 @@ - - - - - - - - - - - - - - -HDAnalyzeR - - - - - - - - - - - - - - - - - - - - - - - - - - -

HDAnalyzeR

- - - -
-

What is HDAnalyzeR?

-

HDAnalyzeR is an R package developed to streamline and enhance -proteomics analysis, particularly for biomarker discovery from blood -plasma samples. It is developed by the Human Disease Blood Atlas group. -This vignette will guide you through the essential steps to use the -package from data loading to biomarker identification. Let’s get started -by loading the package!

-
library(HDAnalyzeR)
-
-
-

How to get Help?

-

If you need help with the package, you can use the ? -operator to access the documentation for a specific function. For -example, to access the documentation for the -hd_initialize() function, you can run:

-
?hd_initialize
-

All package functions can be found in the HDAnalyzeR -Function Reference page. For whole example workflows you can visit -any of the package articles.

-

This document introduces you to HDAnalyzeR’s basic set of tools, and -shows you how to analyze and identify biomarkers in a artificial dataset -of cancer blood plasma samples.

-
-
-

Loading the Data

-

First, we load the package’s example_data and -example_metadata.

-
head(example_data)
-#>      DAid    Sample  OlinkID UniProt  Assay           Panel        NPX
-#> 1 DA00001 AML_syn_1 OID21311  Q9BTE6 AARSD1        Oncology  3.3903461
-#> 2 DA00001 AML_syn_1 OID21280  P00519   ABL1        Oncology  2.7588517
-#> 3 DA00001 AML_syn_1 OID21269  P09110  ACAA1        Oncology  1.7070090
-#> 4 DA00001 AML_syn_1 OID20159  P16112   ACAN Cardiometabolic  0.0332709
-#> 5 DA00001 AML_syn_1 OID20105  Q9BYF1   ACE2 Cardiometabolic  1.7553590
-#> 6 DA00001 AML_syn_1 OID20124  Q15067  ACOX1 Cardiometabolic -0.9192835
-#>   Assay_Warning QC_Warning PlateID
-#> 1          PASS       PASS  Run001
-#> 2          PASS       PASS  Run001
-#> 3          PASS       PASS  Run001
-#> 4          PASS       PASS  Run001
-#> 5          PASS       PASS  Run001
-#> 6          PASS       PASS  Run001
-
head(example_metadata)
-#>      DAid    Sample Disease   Stage Grade Sex Age  BMI Cohort
-#> 1 DA00001 AML_syn_1     AML       2  <NA>   F  42 22.7   UCAN
-#> 2 DA00002 AML_syn_2     AML Unknown  <NA>   M  69 33.1   UCAN
-#> 3 DA00003 AML_syn_3     AML       2  <NA>   F  61 26.2   UCAN
-#> 4 DA00004 AML_syn_4     AML Unknown  <NA>   M  54 28.1   UCAN
-#> 5 DA00005 AML_syn_5     AML       2  <NA>   F  57 21.4   UCAN
-#> 6 DA00006 AML_syn_6     AML Unknown  <NA>   M  86 33.9   UCAN
-

We will initialize the HDAnalyzeR object with the data and metadata. -The hd_initialize() function will create an object that -will be used throughout the analysis. We will specify the sample ID, -variable name, and value name in order to correctly load and widen the -data. In this case, the data is in long format, so we set -is_wide = FALSE.

-
hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-
# See the data transformed into a wide format!
-head(hd_obj$data)
-#> # A tibble: 6 × 101
-#>   DAid    AARSD1  ABL1  ACAA1    ACAN   ACE2  ACOX1   ACP5   ACP6  ACTA2   ACTN4
-#>   <chr>    <dbl> <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
-#> 1 DA00001   3.39  2.76  1.71   0.0333  1.76  -0.919 1.54    2.15   2.81   0.742 
-#> 2 DA00002   1.42  1.25 -0.816 -0.459   0.826 -0.902 0.647   1.30   0.798 -0.0659
-#> 3 DA00003  NA    NA    NA      0.989  NA      0.330 1.37   NA     NA     NA     
-#> 4 DA00004   3.41  3.38  1.69  NA       1.52  NA     0.841   0.582  1.70   0.108 
-#> 5 DA00005   5.01  5.05  0.128  0.401  -0.933 -0.584 0.0265  1.16   2.73   0.350 
-#> 6 DA00006   6.83  1.18 -1.74  -0.156   1.53  -0.721 0.620   0.527  0.772 NA     
-#> # ℹ 90 more variables: ACY1 <dbl>, ADA <dbl>, ADA2 <dbl>, ADAM15 <dbl>,
-#> #   ADAM23 <dbl>, ADAM8 <dbl>, ADAMTS13 <dbl>, ADAMTS15 <dbl>, ADAMTS16 <dbl>,
-#> #   ADAMTS8 <dbl>, ADCYAP1R1 <dbl>, ADGRE2 <dbl>, ADGRE5 <dbl>, ADGRG1 <dbl>,
-#> #   ADGRG2 <dbl>, ADH4 <dbl>, ADM <dbl>, AGER <dbl>, AGR2 <dbl>, AGR3 <dbl>,
-#> #   AGRN <dbl>, AGRP <dbl>, AGXT <dbl>, AHCY <dbl>, AHSP <dbl>, AIF1 <dbl>,
-#> #   AIFM1 <dbl>, AK1 <dbl>, AKR1B1 <dbl>, AKR1C4 <dbl>, AKT1S1 <dbl>,
-#> #   AKT3 <dbl>, ALCAM <dbl>, ALDH1A1 <dbl>, ALDH3A1 <dbl>, ALPP <dbl>, …
-
-

📓 In real-world scenarios, you would load your own data and metadata -files instead of using the example dataset.

-
-
-

⚠️ It is strongly recommended to use hd_initialize() -function to initialize an HDAnalyzeR object before starting any -analysis. This initialization step ensures that your input data and -metadata are correctly formatted and compatible with the package. While -it is possible to use HDAnalyzeR functions on data outside of an -HDAnalyzeR object, users must ensure that the data adheres to the -required structure: the first column must contain sample IDs, followed -by numeric columns in wide format. Additionally, the metadata must -include the same set of sample IDs. Proper initialization or careful -adherence to these requirements is crucial for accurate and efficient -analysis.

-
-
-
-

Quality Control (QC)

-

hd_qc_summary() provides a comprehensive summary of the -input data and metadata.

-

Regarding the data, it will check the column types, calculate the -percentage of NAs in each column and row and plot their distributions, -calculate variable-variable (in this case protein-protein) correlations, -and create a heatmap of these correlations. Users can also specify the -threshold for reporting variable-variable correlations.

-

Regarding the metadata, it will check the column types and calculate -the percentage of NAs in each column and row exactly as before. It will -also create summary visualizations for the metadata variables like “Sex” -and “Age”.

-

Let’s run the QC summary and take a look at some results!

-
qc_data <- hd_qc_summary(hd_obj, 
-                         variable = "Disease", 
-                         palette = list(Disease = "cancers12", Sex = "sex"), 
-                         cor_threshold = 0.7,
-                         verbose = FALSE)
-
-qc_data$data_summary$na_col_hist
-#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
-

-
qc_data$data_summary$cor_heatmap
-

-
qc_data$metadata_summary$Sex
-

-
qc_data$metadata_summary$Age
-#> Picking joint bandwidth of 6.06
-

-
-
-

Imputation and Dimensionality Reduction

-

As we saw from the QC results, the data contains NAs and other issues -that need to be addressed. We will impute missing values using K-nearest -neighbors (KNN) with 5 neighbors via hd_impute_knn().

-
hd_obj <- hd_impute_knn(hd_obj, k = 5, verbose = FALSE)
-head(hd_obj$data)
-#> # A tibble: 6 × 101
-#>   DAid    AARSD1  ABL1  ACAA1    ACAN   ACE2  ACOX1   ACP5  ACP6 ACTA2   ACTN4
-#>   <chr>    <dbl> <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>   <dbl>
-#> 1 DA00001   3.39  2.76  1.71   0.0333  1.76  -0.919 1.54   2.15  2.81   0.742 
-#> 2 DA00002   1.42  1.25 -0.816 -0.459   0.826 -0.902 0.647  1.30  0.798 -0.0659
-#> 3 DA00003   3.70  2.85  2.27   0.989   0.184  0.330 1.37   0.876 1.33   0.570 
-#> 4 DA00004   3.41  3.38  1.69   0.415   1.52   1.23  0.841  0.582 1.70   0.108 
-#> 5 DA00005   5.01  5.05  0.128  0.401  -0.933 -0.584 0.0265 1.16  2.73   0.350 
-#> 6 DA00006   6.83  1.18 -1.74  -0.156   1.53  -0.721 0.620  0.527 0.772  0.446 
-#> # ℹ 90 more variables: ACY1 <dbl>, ADA <dbl>, ADA2 <dbl>, ADAM15 <dbl>,
-#> #   ADAM23 <dbl>, ADAM8 <dbl>, ADAMTS13 <dbl>, ADAMTS15 <dbl>, ADAMTS16 <dbl>,
-#> #   ADAMTS8 <dbl>, ADCYAP1R1 <dbl>, ADGRE2 <dbl>, ADGRE5 <dbl>, ADGRG1 <dbl>,
-#> #   ADGRG2 <dbl>, ADH4 <dbl>, ADM <dbl>, AGER <dbl>, AGR2 <dbl>, AGR3 <dbl>,
-#> #   AGRN <dbl>, AGRP <dbl>, AGXT <dbl>, AHCY <dbl>, AHSP <dbl>, AIF1 <dbl>,
-#> #   AIFM1 <dbl>, AK1 <dbl>, AKR1B1 <dbl>, AKR1C4 <dbl>, AKT1S1 <dbl>,
-#> #   AKT3 <dbl>, ALCAM <dbl>, ALDH1A1 <dbl>, ALDH3A1 <dbl>, ALPP <dbl>, …
-

After imputation, we will run Principal Component Analysis (PCA) via -hd_auto_pca() and Uniform Manifold Approximation and -Projection (UMAP) via hd_auto_umap() to check for outliers, -batch effects, and other potential issues. These are wrapper functions -that will run the necessary steps to perform PCA and UMAP, including -scaling the data and coloring the plots by metadata variables.

-
pca_res <- hd_auto_pca(hd_obj, 
-                       components = 10, 
-                       plot_color = "Disease", 
-                       plot_palette = "cancers12")
-
-head(pca_res$pca_res)
-#> # A tibble: 6 × 11
-#>   DAid     PC1     PC2    PC3    PC4    PC5    PC6    PC7    PC8     PC9    PC10
-#>   <chr>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>   <dbl>
-#> 1 DA00… -3.69  -4.32   -2.30  -3.15  -2.53  -2.78  -2.77   0.603 -3.20   -0.532 
-#> 2 DA00…  4.17  -2.70    2.10  -0.472 -4.47  -1.96  -1.04   1.51  -0.229  -0.0614
-#> 3 DA00… -3.31   4.74    1.41   0.920 -0.588  0.316 -0.103 -0.441  1.36    0.201 
-#> 4 DA00… -4.78   0.447   1.44   0.134 -1.15  -0.284  0.279  3.34   0.0452 -0.866 
-#> 5 DA00… -5.04  -3.69    0.820 -5.73  -0.697 -3.87  -1.01   1.96  -0.970  -1.12  
-#> 6 DA00…  0.429  0.0714 -1.86  -7.80   0.809 -2.89  -0.733  0.381 -0.617   1.32
-pca_res$pca_plot
-#> Ignoring unknown labels:
-#> • Color : "Disease"
-

-
pca_res$pca_loadings_plot
-

-
pca_res$pca_variance_plot
-

-

Let’s color UMAP by Sex instead of Disease to see if there are any -sex related differences.

-
umap_res <- hd_auto_umap(hd_obj, 
-                         plot_color = "Sex", 
-                         plot_palette = "sex")
-
-head(umap_res$umap_res)
-#> # A tibble: 6 × 3
-#>   DAid    UMAP1  UMAP2
-#>   <chr>   <dbl>  <dbl>
-#> 1 DA00001  2.36 -0.770
-#> 2 DA00002 -2.36 -1.34 
-#> 3 DA00003  1.84  1.60 
-#> 4 DA00004  2.44  1.25 
-#> 5 DA00005  2.63 -0.683
-#> 6 DA00006 -1.28  1.22
-umap_res$umap_plot
-#> Ignoring unknown labels:
-#> • Color : "Sex"
-

-
-
-

Biomarker Identification

-
-

Differential Expression Analysis

-

We will run a differential expression analysis to identify potential -biomarkers. We will use hd_de_limma() so that we will be -able to correct also for Sex and Age. This method will help us pinpoint -proteins that are significantly different between conditions. We will -run the analysis 3 times for different case classes against all other -classes each time. We will also run hd_plot_volcano() to -visualize the results. We will present only the results for Acute -Myeloid Leukemia (AML).

-
de_res_aml <- hd_de_limma(hd_obj,
-                          variable = "Disease",
-                          case = "AML",
-                          correct = c("Sex", "Age")) |> 
-  hd_plot_volcano(logfc_lim = 1)
-
-head(de_res_aml$de_res)
-#> # A tibble: 6 × 10
-#>   Feature  logFC   CI.L  CI.R AveExpr     t  P.Value adj.P.Val     B Disease
-#>   <chr>    <dbl>  <dbl> <dbl>   <dbl> <dbl>    <dbl>     <dbl> <dbl> <chr>  
-#> 1 AZU1     1.54   1.23   1.84   0.226  9.97 9.51e-22  9.51e-20  38.5 AML    
-#> 2 ADA      1.40   1.12   1.68   0.924  9.75 6.37e-21  3.19e-19  36.6 AML    
-#> 3 APBB1IP  1.08   0.799  1.35  -0.253  7.62 1.04e-13  3.45e-12  20.3 AML    
-#> 4 APEX1    1.45   1.07   1.82   0.515  7.56 1.50e-13  3.74e-12  19.9 AML    
-#> 5 ANGPT1  -1.63  -2.09  -1.16   1.72  -6.87 1.65e-11  3.29e-10  15.3 AML    
-#> 6 ARTN     0.826  0.572  1.08   0.423  6.38 3.60e-10  6.00e- 9  12.3 AML
-de_res_aml$volcano_plot
-

-
de_res_crc <- hd_de_limma(hd_obj,
-                          variable = "Disease",
-                          case = "CRC",
-                          correct = c("Sex", "Age"))
-
de_res_lungc <- hd_de_limma(hd_obj,
-                            variable = "Disease",
-                            case = "LUNGC",
-                            correct = c("Sex", "Age"))
-

We can also summarize the results via -hd_plot_de_summary(). In order to use this function we need -to store all the results in a list.

-
de_res <- list("AML" = de_res_aml, 
-               "CRC" = de_res_crc, 
-               "LUNGC" = de_res_lungc)
-
-de_summary <- hd_plot_de_summary(de_res, class_palette = "cancers12")
-
de_summary$de_barplot
-

-
de_summary$upset_plot_up
-

-
de_summary$upset_plot_down
-

-
-
-

Machine Learning Classification Model

-

In addition to differential expression analysis, we will use a -regularized machine learning classification model (LASSO) to identify -significant features. This model will help us understand which protein -features are most predictive of the conditions being studied.

-

We will start by creating the data split with -hd_split_data() so we don’t have any data leakage between -train and test data and then we will run the model with -hd_model_rreg(). Once again, we will run the analysis 3 -times but we will present only the results for AML.

-
hd_split <- hd_split_data(hd_obj, variable = "Disease")
-
-lasso_res_aml <- hd_model_rreg(hd_split, 
-                               variable = "Disease",
-                               case = "AML",
-                               mixture = 1,
-                               grid_size = 5,
-                               palette = "cancers12",
-                               verbose = FALSE)
-
-lasso_res_aml$metrics
-#> $accuracy
-#> [1] 0.7959184
-#> 
-#> $sensitivity
-#> [1] 0.9166667
-#> 
-#> $specificity
-#> [1] 0.7851852
-#> 
-#> $auc
-#> [1] 0.9351852
-#> 
-#> $confusion_matrix
-#>           Truth
-#> Prediction   0   1
-#>          0 106   1
-#>          1  29  11
-lasso_res_aml$roc_curve
-

-
lasso_res_aml$probability_plot
-

-
lasso_res_aml$feat_imp_plot
-

-
lasso_res_crc <- hd_model_rreg(hd_split, 
-                               variable = "Disease",
-                               case = "CRC",
-                               mixture = 1,
-                               grid_size = 5,
-                               palette = "cancers12",
-                               verbose = FALSE)
-
lasso_res_lungc <- hd_model_rreg(hd_split, 
-                                 variable = "Disease",
-                                 case = "LUNGC",
-                                 mixture = 1,
-                                 grid_size = 5,
-                                 palette = "cancers12",
-                                 verbose = FALSE)
-

We can get a visual summary of the results via -hd_plot_model_summary(). In order to use this function we -need to store the results in a list as before.

-
lasso_res <- list("AML" = lasso_res_aml, 
-                  "CRC" = lasso_res_crc, 
-                  "LUNGC" = lasso_res_lungc)
-
-features_summary <- hd_plot_model_summary(lasso_res, class_palette = "cancers12")
-
features_summary$metrics_barplot
-#> Ignoring unknown labels:
-#> • colour : "Metric"
-

-
features_summary$features_barplot
-

-
features_summary$upset_plot_features
-

-
-
-
-

One step further

-

The final step involves performing a pathway enrichment analysis to -get more biological information about the features. In this example, we -will run an Over-Representation Analysis for the up-regulated proteins -from differential expression, that are also identified as features by -the ML model. We will use the Gene Ontology (GO) database and show the -results only for AML. We will use hd_ora() and -hd_plot_ora() to perform the analysis and visualize its -results.

-
# Extract the proteins identified by both DE and Lasso
-de_proteins <- de_res_aml$de_res |> 
-  dplyr::filter(adj.P.Val < 0.05) |> 
-  dplyr::pull(Feature)
-
-lasso_proteins <- lasso_res_aml$features |> 
-  dplyr::filter(Scaled_Importance > 0) |> 
-  dplyr::pull(Feature)
-
-intersect_proteins <- intersect(de_proteins, lasso_proteins)
-
-# Perform ORA with GO database and visualize results
-enrichment <- hd_ora(intersect_proteins, 
-                     database = "GO", 
-                     ontology = "BP", 
-                     pval_lim = 0.2)  # We will use a higher p-value threshold for this example
-enrichment_plots <- hd_plot_ora(enrichment)
-
-enrichment_plots$dotplot
-

-
enrichment_plots$treeplot
-

-
enrichment_plots$cnetplot
-

-
-

📓 Remember that these data are a dummy-dataset with artificial data -and the results in this guide should not be interpreted as real results. -The purpose of this vignette is to show you how to use the package and -its functions.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats     graphics  grDevices datasets  utils     methods   base     
-#> 
-#> other attached packages:
-#> [1] HDAnalyzeR_1.0.0
-#> 
-#> loaded via a namespace (and not attached):
-#>   [1] RcppAnnoy_0.0.22        splines_4.5.1           ggplotify_0.1.3        
-#>   [4] tibble_3.3.0            R.oo_1.27.1             hardhat_1.4.2          
-#>   [7] rpart_4.1.24            sparsevctrs_0.3.4       lifecycle_1.0.4        
-#>  [10] globals_0.18.0          lattice_0.22-6          MASS_7.3-65            
-#>  [13] backports_1.5.0         SnowballC_0.7.1         magrittr_2.0.4         
-#>  [16] limma_3.64.3            sass_0.4.10             rmarkdown_2.30         
-#>  [19] jquerylib_0.1.4         yaml_2.3.10             ggtangle_0.0.8         
-#>  [22] askpass_1.2.1           reticulate_1.44.0       cowplot_1.2.0          
-#>  [25] DBI_1.2.3               RColorBrewer_1.1-3      lubridate_1.9.4        
-#>  [28] purrr_1.2.0             R.utils_2.13.0          BiocGenerics_0.54.1    
-#>  [31] yulab.utils_0.2.1       nnet_7.3-20             rappdirs_0.3.3         
-#>  [34] ipred_0.9-15            lava_1.8.2              GenomeInfoDbData_1.2.14
-#>  [37] IRanges_2.42.0          embed_1.2.0             S4Vectors_0.46.0       
-#>  [40] enrichplot_1.28.4       ggrepel_0.9.6           tokenizers_0.3.0       
-#>  [43] listenv_0.10.0          tidytree_0.4.6          pheatmap_1.0.13        
-#>  [46] umap_0.2.10.0           RSpectra_0.16-2         parallelly_1.45.1      
-#>  [49] codetools_0.2-20        DOSE_4.2.0              tidyselect_1.2.1       
-#>  [52] shape_1.4.6.1           aplot_0.2.9             UCSC.utils_1.4.0       
-#>  [55] farver_2.1.2            stats4_4.5.1            jsonlite_2.0.0         
-#>  [58] ggridges_0.5.7          survival_3.8-3          iterators_1.0.14       
-#>  [61] foreach_1.5.2           ggnewscale_0.5.2        tools_4.5.1            
-#>  [64] treeio_1.32.0           Rcpp_1.1.0              glue_1.8.0             
-#>  [67] prodlim_2025.04.28      gridExtra_2.3           xfun_0.54              
-#>  [70] qvalue_2.40.0           tidyheatmaps_0.2.1      GenomeInfoDb_1.44.3    
-#>  [73] dplyr_1.1.4             withr_3.0.2             BiocManager_1.30.25    
-#>  [76] fastmap_1.2.0           openssl_2.3.4           digest_0.6.38          
-#>  [79] timechange_0.3.0        R6_2.6.1                gridGraphics_0.5-1     
-#>  [82] colorspace_2.1-2        GO.db_3.21.0            RSQLite_2.4.4          
-#>  [85] R.methodsS3_1.8.2       UpSetR_1.4.0            utf8_1.2.6             
-#>  [88] tidyr_1.3.1             generics_0.1.4          renv_1.1.4             
-#>  [91] data.table_1.17.8       recipes_1.3.1           class_7.3-23           
-#>  [94] httr_1.4.7              uwot_0.2.4              pkgconfig_2.0.3        
-#>  [97] gtable_0.3.6            parsnip_1.3.3           timeDate_4051.111      
-#> [100] dials_1.4.2             blob_1.2.4              GPfit_1.0-9            
-#> [103] S7_0.2.0                workflows_1.3.0         XVector_0.48.0         
-#> [106] furrr_0.3.1             clusterProfiler_4.16.0  janeaustenr_1.0.0      
-#> [109] htmltools_0.5.8.1       fgsea_1.34.2            scales_1.4.0           
-#> [112] Biobase_2.68.0          png_0.1-8               gower_1.0.2            
-#> [115] ggfun_0.2.0             knitr_1.50              reshape2_1.4.5         
-#> [118] nlme_3.1-168            org.Hs.eg.db_3.21.0     cachem_1.1.0           
-#> [121] stringr_1.6.0           rsample_1.3.1           parallel_4.5.1         
-#> [124] AnnotationDbi_1.70.0    pillar_1.11.1           grid_4.5.1             
-#> [127] vctrs_0.6.5             tune_2.0.1              vip_0.4.1              
-#> [130] lhs_1.2.0               yardstick_1.3.2         evaluate_1.0.5         
-#> [133] cli_3.6.5               compiler_4.5.1          rlang_1.1.6            
-#> [136] crayon_1.5.3            tidytext_0.4.3          future.apply_1.20.0    
-#> [139] labeling_0.4.3          plyr_1.8.9              forcats_1.0.1          
-#> [142] fs_1.6.6                stringi_1.8.7           BiocParallel_1.42.2    
-#> [145] Biostrings_2.76.0       lazyeval_0.2.2          glmnet_4.1-10          
-#> [148] GOSemSim_2.34.0         Matrix_1.7-4            patchwork_1.3.2        
-#> [151] bit64_4.6.0-1           future_1.67.0           ggplot2_4.0.0          
-#> [154] KEGGREST_1.48.1         statmod_1.5.1           tailor_0.1.0           
-#> [157] igraph_2.2.1            broom_1.0.10            memoise_2.0.1          
-#> [160] bslib_0.9.0             ggtree_3.16.3           fastmatch_1.1-6        
-#> [163] bit_4.6.0               gson_0.1.0              ape_5.8-1              
-#> [166] DiceDesign_1.10
-
- - - - - - - - - - - diff --git a/inst/doc/before_using_hdanalyzer.html b/inst/doc/before_using_hdanalyzer.html deleted file mode 100644 index 1429249..0000000 --- a/inst/doc/before_using_hdanalyzer.html +++ /dev/null @@ -1,461 +0,0 @@ - - - - - - - - - - - - - - -Data Preparation for HDAnalyzeR: What You Need Before Using the Package - - - - - - - - - - - - - - - - - - - - - - - - - - -

Data Preparation for HDAnalyzeR: What You -Need Before Using the Package

- - - -
-

Introduction

-

Before you can start using HDAnalyzeR for biomarker discovery, it is -essential to ensure that your data is prepared in the appropriate -format. While HDAnalyzeR offers a variety of powerful tools, it does not -include technology-specific quality control (QC) or preprocessing -functions. This design choice allows the package to be flexible and -usable with a wide range of proteomics technologies without being -limited to specific workflows. Many labs already use their own -preprocessing pipelines, tailored to their specific research needs and -technologies. Integrating all these varied pipelines into the package -would make it unnecessarily complex and restrictive.

-

Therefore, the data provided to HDAnalyzeR must be preprocessed in -order to follow specific requirements (Table 1).

-

Table 1. General requirements

- ---- - - - - - - - - - - - - - - - - - - - - - - - - -
RequirementDescription
DataLong format requires at least three columns: -sample/ID, feature (e.g. proteins, genes, etc.), and measurement -(e.g. TPM, NPX, RFU…).
Wide format requires the -following structure: rows = samples; columns = sample/ID (first column) -+ features.
MetadataTable with one row per sample; includes sample/ID and any required -columns (e.g., group labels, covariates).
Standardized feature IDsRecommended to use gene symbols, UniProt IDs, etc., depending on -platform. Feature IDs may be platform-specific. Mapping tables can help -translate to gene symbols when required (e.g., for pathway -analyses).
Matching sample IDsSample IDs must match exactly between metadata and expression -matrix.
-
-

⚠️ Although HDAnalyzeR is primarily designed for proteomics data, -some functions in the package can also be applied to other omics data -types such as genomics, transcriptomics, or metabolomics. However, it is -important to note that the specific choice of functions and how they -should be used will depend on the study’s goals and the nature of the -data. HDAnalyzeR does not provide explicit guarantees and the user must -make informed decisions regarding its application to other types of -omics data. Documentation and the source code of all functions is freely -available.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats     graphics  grDevices utils     datasets  methods   base     
-#> 
-#> other attached packages:
-#> [1] HDAnalyzeR_1.0.0
-#> 
-#> loaded via a namespace (and not attached):
-#>  [1] jsonlite_2.0.0      dplyr_1.1.4         compiler_4.5.1      BiocManager_1.30.25
-#>  [5] renv_1.1.4          crayon_1.5.3        tidyselect_1.2.1    parallel_4.5.1     
-#>  [9] tidyr_1.3.1         jquerylib_0.1.4     globals_0.18.0      yaml_2.3.10        
-#> [13] fastmap_1.2.0       R6_2.6.1            generics_0.1.4      knitr_1.50         
-#> [17] future_1.67.0       tibble_3.3.0        bslib_0.9.0         pillar_1.11.1      
-#> [21] rlang_1.1.6         utf8_1.2.6          cachem_1.1.0        xfun_0.54          
-#> [25] sass_0.4.10         cli_3.6.5           withr_3.0.2         magrittr_2.0.4     
-#> [29] digest_0.6.38       lifecycle_1.0.4     vctrs_0.6.5         evaluate_1.0.5     
-#> [33] glue_1.8.0          rsample_1.3.1       listenv_0.10.0      furrr_0.3.1        
-#> [37] codetools_0.2-20    parallelly_1.45.1   rmarkdown_2.30      purrr_1.2.0        
-#> [41] tools_4.5.1         pkgconfig_2.0.3     htmltools_0.5.8.1
-
- - - - - - - - - - - diff --git a/inst/doc/classification.html b/inst/doc/classification.html deleted file mode 100644 index 18dacef..0000000 --- a/inst/doc/classification.html +++ /dev/null @@ -1,926 +0,0 @@ - - - - - - - - - - - - - - -Machine Learning Models - - - - - - - - - - - - - - - - - - - - - - - - - - -

Machine Learning Models

- - - -

This vignette will show you how you can easily construct machine -learning pipelines using HDAnalyzeR. We will load HDAnalyzeR and dplyr, -load the example data and metadata that come with the package and -initialize the HDAnalyzeR object.

-
-

Loading the Data

-
library(HDAnalyzeR)
-library(dplyr)
-
-hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-
-

📓 In the whole vignette the verbose parameter of the -model functions will be set to FALSE in order to keep this guide clean -and concise. However, we recommend to leave it to default (TRUE) in -order to know the model’s progress and that everything is running -smoothly.

-
-
-
-

Splitting the Data

-

First, we will create the data split object using the -hd_split_data() function. This function will create a list -of the train and test sets. We can change the ratio of the train and -test sets, the seed for reproducibility, and the metadata variable to -classify. At this stage, we can also add metadata columns as -predictors.

-

We will use the Disease column as the variable to -classify and the Sex and Age columns as a -metadata predictor.

-
split_obj <- hd_split_data(hd_obj, 
-                           variable = "Disease", 
-                           ratio = 0.8, 
-                           seed = 123, 
-                           metadata_cols = c("Sex", "Age"))
-
-
-

Running the Model

-
-

Regularized Regression

-

Let’s start with a regularized regression LASSO model via -hd_model_rreg(). Exactly like in the previous vignette with -the differential expression functions, we have to state the variable, -case and control(s) groups. To do specifically LASSO we will set the -mixture parameter to 1. We will also set the -verbose parameter to FALSE to not print the -progress of the model in shake of clarity for this vignette.

-
model_res <- hd_model_rreg(split_obj,
-                           variable = "Disease",
-                           case = "AML",
-                           control = c("CLL", "MYEL", "GLIOM"),
-                           grid_size = 5,
-                           mixture = 1,
-                           verbose = FALSE)
-
-model_res$final_workflow
-#> ══ Workflow ═══════════════════════════════════════════════════════════════════════
-#> Preprocessor: Recipe
-#> Model: logistic_reg()
-#> 
-#> ── Preprocessor ───────────────────────────────────────────────────────────────────
-#> 5 Recipe Steps
-#> 
-#> • step_dummy()
-#> • step_nzv()
-#> • step_normalize()
-#> • step_corr()
-#> • step_impute_knn()
-#> 
-#> ── Model ──────────────────────────────────────────────────────────────────────────
-#> Logistic Regression Model Specification (classification)
-#> 
-#> Main Arguments:
-#>   penalty = 1.26891735029786e-10
-#>   mixture = mixture
-#> 
-#> Computational engine: glmnet
-model_res$metrics
-#> $accuracy
-#> [1] 0.7352941
-#> 
-#> $sensitivity
-#> [1] 1
-#> 
-#> $specificity
-#> [1] 0.625
-#> 
-#> $auc
-#> [1] 0.9375
-#> 
-#> $confusion_matrix
-#>           Truth
-#> Prediction  0  1
-#>          0 15  0
-#>          1  9 10
-model_res$roc_curve
-

-
model_res$probability_plot
-

-
model_res$feat_imp_plot
-

-

We can change several parameters in the hd_model_rreg() -function. For example, we can change the number of cross-validation -folds, the number of grid points for the hyperparameter optimization, or -the feature correlation threshold. Also, exactly as with the DE -functions, if the control parameter is not set, the -function will use all the other classes as controls. For more -information, please refer to hd_model_rreg() -documentation.

-

We will also set mixture to NULL to allow the model to optimize this -parameter as well (elastic net regression instead of LASSO) and set a -palette for our classes.

-
model_res <- hd_model_rreg(split_obj,
-                           case = "AML",
-                           cv_sets = 3,
-                           grid_size = 5,
-                           cor_threshold = 0.7,
-                           palette = "cancers12",
-                           verbose = FALSE)
-
-model_res$final_workflow
-#> ══ Workflow ═══════════════════════════════════════════════════════════════════════
-#> Preprocessor: Recipe
-#> Model: logistic_reg()
-#> 
-#> ── Preprocessor ───────────────────────────────────────────────────────────────────
-#> 5 Recipe Steps
-#> 
-#> • step_dummy()
-#> • step_nzv()
-#> • step_normalize()
-#> • step_corr()
-#> • step_impute_knn()
-#> 
-#> ── Model ──────────────────────────────────────────────────────────────────────────
-#> Logistic Regression Model Specification (classification)
-#> 
-#> Main Arguments:
-#>   penalty = 0.00355590672132398
-#>   mixture = 0.0638105825171806
-#> 
-#> Computational engine: glmnet
-
-
-

Random Forest

-

We can use a different variable to classify like Sex and -even a different algorithm like random forest via -hd_model_rf(). However, do not forget that we should create -a new split object for this new model. In this case, because the classes -are already balanced, we will set the balance_groups -parameter to FALSE to consider all the samples in the training dataset. -Let’s also remove everything except from number of features and AUC from -the variable importance plot title.

-
split_obj <- hd_split_data(hd_obj, variable = "Sex", ratio = 0.8)
-                               
-model_res <- hd_model_rf(split_obj,
-                       variable = "Sex",
-                       case = "F",
-                       palette = "sex",
-                       cv_sets = 3,
-                       grid_size = 5,
-                       balance_groups = FALSE,
-                       plot_title = c("features", "auc"),
-                       verbose = FALSE)
-
-
-

Logistic Regression

-

If our data have a single predictor, we can use -hd_model_lr() instead of hd_model_rreg() to -perform a logistic regression. Random forest can be used as it was for -multiple predictors.

-
hd_obj_single <- hd_initialize(dat = example_data |> filter(Assay == "ADA"), 
-                               metadata = example_metadata, 
-                               is_wide = FALSE, 
-                               sample_id = "DAid",
-                               var_name = "Assay",
-                               value_name = "NPX")
-
-split_obj <- hd_split_data(hd_obj_single, variable = "Disease", ratio = 0.8)
-
-model_res <- hd_model_lr(split_obj, case = "AML", palette = "cancers12", verbose = FALSE)
-
-
-
-

Visualizing Model Features

-

At this point we should also check how our selected protein features -look in boxplots. We will run a model as before, extract the features, -select the top-9 of them based on their importance in the model and plot -them with hd_plot_feature_boxplot(). We can either plot -case vs control or case vs all other classes by changing the -type argument.

-
-

⚠️ In case you have metadata variables as features, you will have to -remove them from the feature vector before using the -hd_plot_feature_boxplot() function as it is made to -visualize protein features.

-
-
hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-
-split_obj <- hd_split_data(hd_obj, variable = "Disease", ratio = 0.8)
-
-model_res <- hd_model_rreg(split_obj, case = "AML", cv_sets = 3, grid_size = 5, verbose = FALSE)
-
-features <- model_res$features |> arrange(desc(Scaled_Importance)) |> head(9) |> pull(Feature)
-
-hd_plot_feature_boxplot(hd_obj, 
-                        features = features, 
-                        case = "AML", 
-                        palette = "cancers12", 
-                        type = "case_vs_control",
-                        points = FALSE)
-

-

-hd_plot_feature_boxplot(hd_obj, 
-                        features = features, 
-                        case = "AML", 
-                        palette = "cancers12", 
-                        type = "case_vs_all")
-

-
-
-

Multi-classification Model

-

We can also do multiclassification predictions with all available -classes in the data. The only thing that we should change is set the -case argument to NULL so that the model understands that we -want to classify all the classes. Let’s see an example with regularized -regression!

-
model_res <- hd_model_rreg(split_obj, 
-                           case = NULL, 
-                           cv_sets = 3, 
-                           grid_size = 5, 
-                           palette = "cancers12",
-                           verbose = FALSE)
-
-model_res$final_workflow
-#> ══ Workflow ═══════════════════════════════════════════════════════════════════════
-#> Preprocessor: Recipe
-#> Model: multinom_reg()
-#> 
-#> ── Preprocessor ───────────────────────────────────────────────────────────────────
-#> 5 Recipe Steps
-#> 
-#> • step_dummy()
-#> • step_nzv()
-#> • step_normalize()
-#> • step_corr()
-#> • step_impute_knn()
-#> 
-#> ── Model ──────────────────────────────────────────────────────────────────────────
-#> Multinomial Regression Model Specification (classification)
-#> 
-#> Main Arguments:
-#>   penalty = 0.00445526787557566
-#>   mixture = 0.784269891628064
-#> 
-#> Computational engine: glmnet
-model_res$roc_curve
-

-
model_res$probability_plot
-

-
model_res$feat_imp_plot
-

-
-
-

Regression instead of Classification

-

Instead of a classification we can run a regression model. That means -that we will try to predict a continuous variable instead of a -categorical one. We can use either hd_model_rreg() or -hd_model_rf() functions with the case -parameter set to NULL. Let’s see an example with the Age -variable. Do not forget that we have to create a new split object for -this new model with Age as the variable of interest.

-
-

⚠️ We should not forget to update the plot_title -argument by changing the metrics from “accuracy”, “sensitivity”, -“apwcificity”, and “auc” to “rmse” and “rsq”.

-
-
split_obj <- hd_split_data(hd_obj, variable = "Age", ratio = 0.8)
-
-model_res <- hd_model_rreg(split_obj, 
-                           variable = "Age",
-                           case = NULL, 
-                           cv_sets = 3, 
-                           grid_size = 2,
-                           plot_title = c("rmse", "rsq", "features", "mixture"),
-                           verbose = FALSE)
-
-model_res$final_workflow
-#> ══ Workflow ═══════════════════════════════════════════════════════════════════════
-#> Preprocessor: Recipe
-#> Model: linear_reg()
-#> 
-#> ── Preprocessor ───────────────────────────────────────────────────────────────────
-#> 5 Recipe Steps
-#> 
-#> • step_dummy()
-#> • step_nzv()
-#> • step_normalize()
-#> • step_corr()
-#> • step_impute_knn()
-#> 
-#> ── Model ──────────────────────────────────────────────────────────────────────────
-#> Linear Regression Model Specification (regression)
-#> 
-#> Main Arguments:
-#>   penalty = 4.45590449619826e-06
-#>   mixture = 0.220184963848442
-#> 
-#> Computational engine: glmnet
-model_res$comparison_plot
-

-
model_res$feat_imp_plot
-

-
-
-

Test the Model on new Data

-

Furthermore, we can validate our trained model in new data. For this -example we will not use another dataset, but we will split the data -initially to create a train and a validation set and then split the -train set to an inner train and a test set. We will use this second -split to initially train the model and then evaluate it with the -validation data. In a real case scenario, you can do either this, or use -a completely different dataset to check that the model generalizes -properly. We will use the hd_model_test() function to do -this. Let’s see an example with the AML model.

-
# Split the data for training and validation sets
-dat <- hd_obj$data
-train_indices <- sample(1:nrow(dat), size = floor(0.8 * nrow(dat)))
-train_data <- dat[train_indices, ]
-validation_data <- dat[-train_indices, ]
-
-hd_object_train <- hd_initialize(train_data, example_metadata, is_wide = TRUE)
-hd_object_val <- hd_initialize(validation_data, example_metadata, is_wide = TRUE)
-
-# Split the training set into training and inner test sets
-split_obj <- hd_split_data(hd_object_train, variable = "Disease")
-
-# Run the regularized regression model pipeline
-model_object <- hd_model_rreg(split_obj,
-                              variable = "Disease",
-                              case = "AML",
-                              grid_size = 2,
-                              palette = "cancers12")
-
-# Run the model evaluation pipeline
-model_res <- hd_model_test(model_object, 
-                           hd_object_train, 
-                           hd_object_val, 
-                           case = "AML", 
-                           palette = "cancers12")
-
-model_res$metrics
-#> $accuracy
-#> [1] 0.8205128
-#> 
-#> $sensitivity
-#> [1] 1
-#> 
-#> $specificity
-#> [1] 0.8018868
-#> 
-#> $auc
-#> [1] 0.9502573
-#> 
-#> $confusion_matrix
-#>           Truth
-#> Prediction  0  1
-#>          0 85  0
-#>          1 21 11
-model_res$test_metrics  # Results from the validation set
-#> $accuracy
-#> [1] 0.8050847
-#> 
-#> $sensitivity
-#> [1] 0.7777778
-#> 
-#> $specificity
-#> [1] 0.8073394
-#> 
-#> $auc
-#> [1] 0.8674822
-#> 
-#> $confusion_matrix
-#>           Truth
-#> Prediction  0  1
-#>          0 88  2
-#>          1 21  7
-model_res$roc_curve
-

-
model_res$test_roc_curve  # Results from the validation set
-

-
-
-

Summarizing Results from Multiple Binary Models

-

To summarize the results for multiple binary models we can use the -hd_plot_model_summary() function. We can create models of -different cases and compare them. Let’s run three different models for -three different cancers and summarize them.

-
-

📓 Do not forget that Ovarian Cancer is sex specific and we should -consider run the analysis only with samples of that sex. We can easily -integrate that into our pipeline using the hd_filter() -function.

-
-
split_obj <- hd_split_data(hd_obj, variable = "Disease")
-
-model_aml <- hd_model_rreg(split_obj, case = "AML", cv_sets = 3, grid_size = 5, verbose = FALSE)
-
-model_gliom <- hd_model_rreg(split_obj, case = "GLIOM", cv_sets = 3, grid_size = 5, verbose = FALSE)
-
-split_obj_sex <- hd_split_data(hd_obj |> hd_filter(variable = "Sex", values = "F", flag = "k"),
-                               variable = "Disease",
-                               ratio = 0.8)
-
-model_ovc <- hd_model_rreg(split_obj_sex, case = "OVC", cv_sets = 3, grid_size = 5, verbose = FALSE)
-
model_summary_res <- hd_plot_model_summary(list("AML" = model_aml, 
-                                                "GLIOM" = model_gliom, 
-                                                "OVC" = model_ovc), 
-                                           class_palette = "cancers12")
-
model_summary_res$metrics_barplot
-#> Ignoring unknown labels:
-#> • colour : "Metric"
-

-
model_summary_res$features_barplot
-

-
model_summary_res$upset_plot_features
-

-

In case we have one case and multiple controls we can use the -hd_plot_feature_heatmap() function to visualize the protein -features in a heatmap. This function is useful as we can easily see if -the same features are important in multiple models. Let’s see an example -with the AML model and 3 different controls groups. We will combine DE -results of the same comparisons.

-
model_cll <- hd_model_rreg(split_obj, case = "AML", control = "CLL", cv_sets = 3, grid_size = 5, verbose = FALSE)
-
-model_blood <- hd_model_rreg(split_obj, 
-                             case = "AML", 
-                             control = c("CLL", "MYEL", "LYMPH"), 
-                             cv_sets = 3, 
-                             grid_size = 5, 
-                             verbose = FALSE)
-
-model_all <- hd_model_rreg(split_obj, case = "AML", cv_sets = 3, grid_size = 5, verbose = FALSE)
-
-de_cll <- hd_de_limma(hd_obj, case = "AML", control = "CLL", correct = c("Sex", "Age"))
-
-de_blood <- hd_de_limma(hd_obj, 
-                              case = "AML", 
-                              control = c("CLL", "MYEL", "LYMPH"), 
-                              correct = c("Sex", "Age"))
-
-de_all <- hd_de_limma(hd_obj, case = "AML", correct = c("Sex", "Age"))
-
hd_plot_feature_heatmap(de_results = list("CLL" = de_cll, 
-                                          "Blood" = de_blood, 
-                                          "All" = de_all), 
-                        model_results = list("CLL" = model_cll, 
-                                             "Blood" = model_blood, 
-                                             "All" = model_all), 
-                        order_by = "CLL")
-

-

Finally, we can use the hd_plot_feature_network() -function to visualize the protein features in a network. This function -is useful as we can easily see the connections between the features and -the importance of each feature in the model. Let’s see an example with -the same 3 models from before.

-
feature_panel <- model_aml[["features"]] |>
-  filter(Scaled_Importance > 0.5) |>
-  mutate(Class = "AML") |>
-  bind_rows(model_gliom[["features"]] |>
-              filter(Scaled_Importance > 0.5) |>
-              mutate(Class = "GLIOM"),
-            model_ovc[["features"]] |>
-              filter(Scaled_Importance > 0.5) |>
-              mutate(Class = "OVC"))
-
-print(head(feature_panel))  # Preview of the feature panel
-#> # A tibble: 6 × 5
-#>   Feature  Importance Sign  Scaled_Importance Class
-#>   <fct>         <dbl> <chr>             <dbl> <chr>
-#> 1 ANGPT1        1.31  NEG               1     AML  
-#> 2 ADGRG1        1.02  POS               0.779 AML  
-#> 3 AMY2A         0.848 POS               0.648 AML  
-#> 4 ADAMTS16      0.788 NEG               0.602 AML  
-#> 5 ADA           0.769 POS               0.587 AML  
-#> 6 ADAM8         0.721 NEG               0.551 AML
-
-hd_plot_feature_network(feature_panel,
-                        plot_color = "Scaled_Importance",
-                        class_palette = "cancers12")
-

-
-

📓 Remember that these data are a dummy-dataset with artificial data -and the results in this guide should not be interpreted as real results. -The purpose of this vignette is to show you how to use the package and -its functions.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats     graphics  grDevices datasets  utils     methods   base     
-#> 
-#> other attached packages:
-#> [1] dplyr_1.1.4      HDAnalyzeR_1.0.0
-#> 
-#> loaded via a namespace (and not attached):
-#>   [1] gridExtra_2.3       rlang_1.1.6         magrittr_2.0.4     
-#>   [4] furrr_0.3.1         tailor_0.1.0        compiler_4.5.1     
-#>   [7] vctrs_0.6.5         lhs_1.2.0           stringr_1.6.0      
-#>  [10] tune_2.0.1          pkgconfig_2.0.3     shape_1.4.6.1      
-#>  [13] fastmap_1.2.0       ggraph_2.2.2        labeling_0.4.3     
-#>  [16] utf8_1.2.6          rmarkdown_2.30      prodlim_2025.04.28 
-#>  [19] ggbeeswarm_0.7.2    UpSetR_1.4.0        purrr_1.2.0        
-#>  [22] xfun_0.54           glmnet_4.1-10       cachem_1.1.0       
-#>  [25] jsonlite_2.0.0      recipes_1.3.1       tweenr_2.0.3       
-#>  [28] vip_0.4.1           parallel_4.5.1      R6_2.6.1           
-#>  [31] bslib_0.9.0         rsample_1.3.1       stringi_1.8.7      
-#>  [34] RColorBrewer_1.1-3  ranger_0.17.0       limma_3.64.3       
-#>  [37] parallelly_1.45.1   boot_1.3-32         rpart_4.1.24       
-#>  [40] lubridate_1.9.4     jquerylib_0.1.4     Rcpp_1.1.0         
-#>  [43] dials_1.4.2         iterators_1.0.14    knitr_1.50         
-#>  [46] future.apply_1.20.0 zoo_1.8-14          multiROC_1.1.1     
-#>  [49] igraph_2.2.1        Matrix_1.7-4        splines_4.5.1      
-#>  [52] nnet_7.3-20         timechange_0.3.0    tidyselect_1.2.1   
-#>  [55] viridis_0.6.5       yaml_2.3.10         timeDate_4051.111  
-#>  [58] codetools_0.2-20    listenv_0.10.0      lattice_0.22-6     
-#>  [61] tibble_3.3.0        plyr_1.8.9          withr_3.0.2        
-#>  [64] S7_0.2.0            evaluate_1.0.5      future_1.67.0      
-#>  [67] survival_3.8-3      polyclip_1.10-7     pillar_1.11.1      
-#>  [70] BiocManager_1.30.25 renv_1.1.4          foreach_1.5.2      
-#>  [73] generics_0.1.4      ggplot2_4.0.0       scales_1.4.0       
-#>  [76] globals_0.18.0      class_7.3-23        glue_1.8.0         
-#>  [79] tools_4.5.1         ggnewscale_0.5.2    data.table_1.17.8  
-#>  [82] gower_1.0.2         forcats_1.0.1       graphlayouts_1.2.2 
-#>  [85] tidygraph_1.3.1     grid_4.5.1          yardstick_1.3.2    
-#>  [88] tidyr_1.3.1         ipred_0.9-15        ggforce_0.5.0      
-#>  [91] beeswarm_0.4.0      vipor_0.4.7         cli_3.6.5          
-#>  [94] DiceDesign_1.10     workflows_1.3.0     parsnip_1.3.3      
-#>  [97] viridisLite_0.4.2   lava_1.8.2          gtable_0.3.6       
-#> [100] GPfit_1.0-9         sass_0.4.10         digest_0.6.38      
-#> [103] ggrepel_0.9.6       farver_2.1.2        memoise_2.0.1      
-#> [106] htmltools_0.5.8.1   lifecycle_1.0.4     hardhat_1.4.2      
-#> [109] statmod_1.5.1       MASS_7.3-65         sparsevctrs_0.3.4
-
- - - - - - - - - - - diff --git a/inst/doc/data_imputation.html b/inst/doc/data_imputation.html deleted file mode 100644 index 80ea3b0..0000000 --- a/inst/doc/data_imputation.html +++ /dev/null @@ -1,628 +0,0 @@ - - - - - - - - - - - - - - -Data Imputation - - - - - - - - - - - - - - - - - - - - - - - - - - -

Data Imputation

- - - -

This vignette will guide you to the different imputation methods -HDAnalyzeR offers. First of all, we will load the package, as well as -dplyr, ggplot2 and patchwork for data manipulation and -visualization.

-
library(HDAnalyzeR)
-library(dplyr)
-library(ggplot2)
-library(patchwork)
-
-

Loading the Data

-

Let’s start with loading the example data and metadata that come with -the package and initialize the HDAnalyzeR object.

-
hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-
-
-

Explore Missing Values

-

We can simply check our data for NA values by using the -hd_qc_summary() as we did in previous vignettes. This time -we will use something specific to NA values, the -hd_na_search() function. This function will return a -summary heatmap showing the distribution of NA values across the data -and metadata variables. This function is ideal to dive into the missing -values and understand if there are any patterns in the missing data. -This is important in order to decide how to handle them (e.g., impute or -remove).

-
na_res <- hd_na_search(hd_obj,
-                       annotation_vars = c("Sex", "Age", "Disease"),
-                       palette = list(Disease = "cancers12",
-                                      Sex = "sex"),
-                       x_labels = FALSE,
-                       y_labels = FALSE)
-
-na_res$na_heatmap
-

-

In this case, we can see that the NA values are generally spread -across the different Assays, samples and metadata variables. There is a -higher concentration of missing values in Myeloma that may require -further investigation. In our case, we will try impute them!

-
-
-

Imputation Methods

-
-

Median Imputation

-

We will start the imputation with the simplest and fastest method, -which is the median imputation by using the -hd_impute_median(). After the imputation, we will check the -sample distribution of a random Assay that contains missing values to -see if these values are imputed logically. In a real case, this check -should be done to more than just one assay.

-
imputed_hd_obj <- hd_impute_median(hd_obj, verbose = FALSE)
-
-plot_before <- hd_obj$data |> 
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightblue", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "Before imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_after <- imputed_hd_obj$data |> 
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightpink", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "After median imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_before + plot_after
-

-

As observed in the plots, the distribution of the ACTA2 assay shifts -after imputation, with an exaggerated median value in the imputed data. -This highlights a key drawback of median imputation: the more missing -values there are, the greater the potential bias.

-
-
-

KNN Imputation

-

A better approach is to use the hd_impute_knn() with 5 -neighbors, which imputes missing values based on the 5-nearest -neighbors. We will use the same assay to compare the imputed data with -the original data.

-
imputed_hd_obj <- hd_impute_knn(hd_obj, k = 5, verbose = FALSE)
-
-plot_before <- hd_obj$data |> 
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightblue", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "Before imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_after <- imputed_hd_obj$data |>
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightpink", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "After KNN imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_before + plot_after
-

-

In this case, the distribution of the ACTA2 assay after imputation is -more similar to the original distribution. This is because the KNN -imputation method uses the nearest neighbors to impute missing values, -which is more accurate and representative than median imputation.

-
-
-

MissForest Imputation

-

Finally, we will use the hd_impute_missForest() method, -which uses the random forest algorithm to impute missing values. We will -use the default values for the number of trees and the number of -iterations.

-
imputed_hd_obj <- hd_impute_missForest(hd_obj, verbose = FALSE)
-
-plot_before <- hd_obj$data |> 
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightblue", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "Before imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_after <- imputed_hd_obj$data |>
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightpink", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "After MissForest imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_before + plot_after
-

-

The MissForest imputation method is usually the most accurate and -also very robust, as it uses the complex random forest algorithm to -impute missing values. This method is particularly useful for large -datasets with complex relationships between variables. On the other -hand, it is by far the most computationally expensive and it would help -parallelize it. You can do that by creating and registering a cluster -with a package like doParallel and then setting the -parallelize argument to “forests” or “variables”.

-
-

📓 All methods assume that the data is missing at random, which is a -common assumption in imputation methods. If the data are missing in a -biased way (either technical or biological), the imputation methods may -introduce bias into the data. In such cases, it is important to -carefully consider the way the data were collected and what they -represent.

-
-
-
-
-

Removing Missing Values instead of Imputing

-

If for any reason you do not want to impute the data, you can use the -hd_omit_na() function to easily remove the rows with -missing values in specific variables. In this example, we will remove -all rows with missing values in any of the assays.

-
imputed_hd_obj <- hd_omit_na(hd_obj)
-
-plot_before <- hd_obj$data |> 
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightblue", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "Before imputation",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_after <- imputed_hd_obj$data |>
-  ggplot(aes(x = ACTA2)) +
-  geom_histogram(fill = "lightpink", color = "black", alpha = 0.5, bins = 30) +
-  labs(title = "After removing missing values",
-       x = "ACTA2", y = "Count") +
-  ylim(0, 90) +
-  theme_minimal()
-
-plot_before + plot_after
-

-

-# Data after removing missing values only in specific columns
-res <- hd_omit_na(hd_obj, columns = "AARSD1")
-res$data
-#> # A tibble: 552 × 101
-#>    DAid    AARSD1  ABL1  ACAA1    ACAN    ACE2  ACOX1    ACP5    ACP6 ACTA2
-#>    <chr>    <dbl> <dbl>  <dbl>   <dbl>   <dbl>  <dbl>   <dbl>   <dbl> <dbl>
-#>  1 DA00001   3.39 2.76   1.71   0.0333  1.76   -0.919  1.54    2.15   2.81 
-#>  2 DA00002   1.42 1.25  -0.816 -0.459   0.826  -0.902  0.647   1.30   0.798
-#>  3 DA00004   3.41 3.38   1.69  NA       1.52   NA      0.841   0.582  1.70 
-#>  4 DA00005   5.01 5.05   0.128  0.401  -0.933  -0.584  0.0265  1.16   2.73 
-#>  5 DA00006   6.83 1.18  -1.74  -0.156   1.53   -0.721  0.620   0.527  0.772
-#>  6 DA00008   2.78 0.812 -0.552  0.982  -0.101  -0.304  0.376  -0.826  1.52 
-#>  7 DA00009   4.39 3.34  -0.452 -0.868   0.395   1.71   1.49   -0.0285 0.200
-#>  8 DA00010   1.83 1.21  -0.912 -1.04   -0.0918 -0.304  1.69    0.0920 2.04 
-#>  9 DA00011   3.48 4.96   3.50  -0.338   4.48    1.26   2.18    1.62   1.79 
-#> 10 DA00012   4.31 0.710 -1.44  -0.218  -0.469  -0.361 -0.0714 -1.30   2.86 
-#> # ℹ 542 more rows
-#> # ℹ 91 more variables: ACTN4 <dbl>, ACY1 <dbl>, ADA <dbl>, ADA2 <dbl>,
-#> #   ADAM15 <dbl>, ADAM23 <dbl>, ADAM8 <dbl>, ADAMTS13 <dbl>, ADAMTS15 <dbl>,
-#> #   ADAMTS16 <dbl>, ADAMTS8 <dbl>, ADCYAP1R1 <dbl>, ADGRE2 <dbl>, ADGRE5 <dbl>,
-#> #   ADGRG1 <dbl>, ADGRG2 <dbl>, ADH4 <dbl>, ADM <dbl>, AGER <dbl>, AGR2 <dbl>,
-#> #   AGR3 <dbl>, AGRN <dbl>, AGRP <dbl>, AGXT <dbl>, AHCY <dbl>, AHSP <dbl>,
-#> #   AIF1 <dbl>, AIFM1 <dbl>, AK1 <dbl>, AKR1B1 <dbl>, AKR1C4 <dbl>, …
-

In this vignette we showed that via HDAnalyzeR you can impute your -data with different methods, each of them with its own advantages and -drawbacks. You can choose the method that best fits your data and your -analysis needs. When using KNN or MissForest imputation methods, you -should experiment with the parameters and look at the distributions of -assays before and after to pick the most suitable.

-
-

📓 Remember that these data are a dummy-dataset with artificial data -and the results in this guide should not be interpreted as real results. -The purpose of this vignette is to show you how to use the package and -its functions.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats     graphics  grDevices datasets  utils     methods   base     
-#> 
-#> other attached packages:
-#> [1] patchwork_1.3.2  ggplot2_4.0.0    dplyr_1.1.4      HDAnalyzeR_1.0.0
-#> 
-#> loaded via a namespace (and not attached):
-#>  [1] tidyselect_1.2.1     timeDate_4051.111    farver_2.1.2        
-#>  [4] S7_0.2.0             fastmap_1.2.0        digest_0.6.38       
-#>  [7] rpart_4.1.24         timechange_0.3.0     lifecycle_1.0.4     
-#> [10] survival_3.8-3       magrittr_2.0.4       compiler_4.5.1      
-#> [13] rlang_1.1.6          sass_0.4.10          rngtools_1.5.2      
-#> [16] tools_4.5.1          utf8_1.2.6           yaml_2.3.10         
-#> [19] data.table_1.17.8    knitr_1.50           labeling_0.4.3      
-#> [22] doRNG_1.8.6.2        RColorBrewer_1.1-3   withr_3.0.2         
-#> [25] purrr_1.2.0          itertools_0.1-3      nnet_7.3-20         
-#> [28] grid_4.5.1           sparsevctrs_0.3.4    future_1.67.0       
-#> [31] globals_0.18.0       scales_1.4.0         iterators_1.0.14    
-#> [34] MASS_7.3-65          cli_3.6.5            rmarkdown_2.30      
-#> [37] generics_0.1.4       future.apply_1.20.0  cachem_1.1.0        
-#> [40] splines_4.5.1        parallel_4.5.1       ggplotify_0.1.3     
-#> [43] BiocManager_1.30.25  yulab.utils_0.2.1    vctrs_0.6.5         
-#> [46] hardhat_1.4.2        Matrix_1.7-4         jsonlite_2.0.0      
-#> [49] gridGraphics_0.5-1   listenv_0.10.0       foreach_1.5.2       
-#> [52] gower_1.0.2          tidyr_1.3.1          jquerylib_0.1.4     
-#> [55] recipes_1.3.1        missForest_1.6.1     glue_1.8.0          
-#> [58] parallelly_1.45.1    codetools_0.2-20     lubridate_1.9.4     
-#> [61] gtable_0.3.6         tibble_3.3.0         pillar_1.11.1       
-#> [64] rappdirs_0.3.3       htmltools_0.5.8.1    ipred_0.9-15        
-#> [67] randomForest_4.7-1.2 lava_1.8.2           R6_2.6.1            
-#> [70] Rdpack_2.6.4         evaluate_1.0.5       lattice_0.22-6      
-#> [73] rbibutils_2.4        pheatmap_1.0.13      renv_1.1.4          
-#> [76] bslib_0.9.0          class_7.3-23         Rcpp_1.1.0          
-#> [79] prodlim_2025.04.28   ranger_0.17.0        xfun_0.54           
-#> [82] fs_1.6.6             tidyheatmaps_0.2.1   pkgconfig_2.0.3
-
- - - - - - - - - - - diff --git a/inst/doc/data_preprocessing.html b/inst/doc/data_preprocessing.html deleted file mode 100644 index aace704..0000000 --- a/inst/doc/data_preprocessing.html +++ /dev/null @@ -1,590 +0,0 @@ - - - - - - - - - - - - - - -Exploratory Data Analysis - - - - - - - - - - - - - - - - - - - - - - - - - - -

Exploratory Data Analysis

- - - -

This vignette will guide you through the initial steps of a -proteomics pipeline, which is the exploratory data analysis. Except from -HDAnalyzeR we will import viridis for the color palettes.

-
library(HDAnalyzeR)
-library(viridis)
-
-

📓 This vignette is just a basic example of how to explore the data -and perform dimensionality reduction to identify potential issues and -trends. However, every dataset is different and may require different -exploration techniques other than the ones contained in the package.

-
-
-

⚠️ Please make sure to preprocess the data and metadata before -proceeding with the analysis. For the most common proteomics data -formats, we provide some recommendations in the package documentation -under articles “Data Preparation for HDAnalyzeR: What You Need Before -Using the Package”.

-
-
-

Loading the Data

-

Let’s start with loading the example data and metadata that come with -the package and initialize the HDAnalyzeR object.

-
hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-
-
-

Exploratory Data Analysis

-

Now, we will perform an automated exploratory data analysis (QC -check) on the data. hd_qc_summary() will return the -percentage of missing values for each column and row and histograms of -their distributions both for data and metadata, the variable-variable -(in this case protein-protein) correlation and the metadata variable -distributions.

-
-

⚠️ This should not be confused with the quality control of the data, -which is a different process that is usually performed right after the -proteomics experiments. The exploratory data analysis is a first step to -understand the data and identify potential issues.

-
-

As variable we should pass the column name of the -metadata that contains the different classes, for example the case and -control groups, the different diseases, or the different stages of a -disease. The palettes are optional and can be used to color the plots of -the metadata variable distributions. You should pass a list with the -variable name and the palette name, for example list(Sex = c(“F” = -“red”, “M” = “blue”), Disease = “cancers12”). As you can see, the -palettes can be either a named vector or a character vector with the -name of a palette from the package. You can display all available -HDAnalyzeR palettes by running hd_show_palettes().

-
qc_data <- hd_qc_summary(hd_obj, 
-                         variable = "Disease", 
-                         palette = list(Disease = "cancers12", Sex = "sex"), 
-                         cor_threshold = 0.7)
-
-qc_data$data_summary$na_col_hist
-

-
qc_data$data_summary$na_row_hist
-

-
qc_data$data_summary$cor_heatmap
-

-
qc_data$metadata_summary$na_col_hist
-

-
qc_data$metadata_summary$na_row_hist
-

-
qc_data$metadata_summary$Sex
-

-
qc_data$metadata_summary$Stage
-

-
qc_data$metadata_summary$Age
-

-
qc_data$metadata_summary$BMI
-

-

From the EDA results we can see that all assays have less than 10% -missing values, while only a few samples have more than 20% missing -values. The protein-protein correlation plot shows that there are no -extremely highly correlated proteins. Also, the metadata variable -Stage has a high percentage of missing values that are not -correctly stated as NA values but as “” or “Unknown”. These are only a -few examples of the information that we can extract from the EDA -results.

-

Based on these we can take decisions on filtering samples with high -missing values, excluding whole assays and metadata variables from -statistical analysis.

-
-
-

Dimensionality Reduction

-

Finally we will run a PCA and UMAP analysis on the data to check if -the data contain any outliers or cluster in an unexpected way. From the -PCA analysis we can also check how the variance in the data is explained -by the different principal components. In this example, we will color -the PCA plot based on Disease, while the UMAP plot based on -Age.

-

The hd_pca() and hd_umap() functions will -run the respective analysis while the hd_plot_dim() will -plot the data on a 2D plane. For the PCA analysis we can also use -hd_plot_pca_loadings() and -hd_plot_pca_variance() to plot the loadings of the -different proteins on the principal components and the variance -explained by each principal component respectively.

-

In the hd_plot_dim() function we can pass the -color parameter to color the data points based on a -metadata variable. We should not forget to pass again the hd_obj object -to the metadata argument so the package can search the -metadata and identify the color variable.

-
pca_res <- hd_pca(hd_obj, components = 15) |> 
-  hd_plot_dim(hd_obj, "PC1", "PC2", color = "Disease", palette = "cancers12", axis_variance = TRUE) |> 
-  hd_plot_pca_loadings(displayed_pcs = 6, displayed_features = 10) |> 
-  hd_plot_pca_variance()
-
-head(pca_res$pca_res)
-#> # A tibble: 6 × 16
-#>   DAid     PC1     PC2    PC3    PC4    PC5    PC6     PC7    PC8    PC9    PC10
-#>   <chr>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
-#> 1 DA00… -3.67  -4.28   -2.34  -3.10  -2.65  -2.73  -2.78    0.433 -3.23  -0.399 
-#> 2 DA00…  4.11  -2.64    2.04  -0.441 -4.43  -1.91  -0.897   1.54  -0.283  0.166 
-#> 3 DA00… -3.34   4.72    1.41   0.881 -0.561  0.308 -0.0612 -0.267  1.42   0.0840
-#> 4 DA00… -4.78   0.443   1.41   0.107 -1.10  -0.262  0.350   3.33  -0.336 -0.748 
-#> 5 DA00… -4.98  -3.67    0.711 -5.70  -0.807 -3.77  -0.969   1.83  -1.29  -0.991 
-#> 6 DA00…  0.395  0.0572 -1.90  -7.75   0.707 -2.70  -0.681   0.301 -0.486  1.32  
-#> # ℹ 5 more variables: PC11 <dbl>, PC12 <dbl>, PC13 <dbl>, PC14 <dbl>,
-#> #   PC15 <dbl>
-pca_res$pca_plot
-

-
pca_res$pca_loadings_plot
-

-
pca_res$pca_variance_plot
-

-

This time we want to color based on a continuous variable, the -Age. In this case, we can either bin the data into -categories using hd_bin_columns() and use a categorical -palette as before or use a continuous color palette directly after the -plot is created. We will do the second by setting the palette argument -to NULL and using the viridis palette afterwards.

-
umap_res <- hd_umap(hd_obj, components = 2) |> 
-  hd_plot_dim(hd_obj, "UMAP1", "UMAP2", color = "Age", palette = NULL)
-
-umap_res$umap_plot + scale_color_viridis()
-

-

We can also run a UMAP where the data points are Proteins (or any -feature) instead of samples. The only thing we need to do is to set the -by_sample parameter to FALSE and pass the -Assay column as the plot_color parameter. We -can also pass a custom palette to color the different assays. This time -we will use the wrapper function hd_auto_umap(). The -wrapper functions are a bit less flexible but they are easier to use and -require less code.

-
umap_res <- hd_auto_umap(hd_obj, 
-                         by_sample = FALSE, 
-                         plot_color = "Assay",
-                         plot_palette = c("ADA" = "darkblue", 
-                                          "ABL1" = "red3", 
-                                          "ACAN" = "green3"))
-
-head(umap_res$umap_res)
-#> # A tibble: 6 × 3
-#>   Assay   UMAP1  UMAP2
-#>   <chr>   <dbl>  <dbl>
-#> 1 AARSD1  3.11  -2.55 
-#> 2 ABL1    1.20  -2.99 
-#> 3 ACAA1  -0.113 -2.16 
-#> 4 ACAN   -1.55   0.184
-#> 5 ACE2   -0.489 -0.668
-#> 6 ACOX1  -0.283  2.51
-umap_res$umap_plot
-

-
-

📓 Remember that these data are a dummy-dataset with artificial data -and the results in this guide should not be interpreted as real results. -The purpose of this vignette is to show you how to use the package and -its functions.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats     graphics  grDevices datasets  utils     methods   base     
-#> 
-#> other attached packages:
-#> [1] viridis_0.6.5     viridisLite_0.4.2 HDAnalyzeR_1.0.0 
-#> 
-#> loaded via a namespace (and not attached):
-#>  [1] tidyselect_1.2.1    timeDate_4051.111   dplyr_1.1.4        
-#>  [4] farver_2.1.2        S7_0.2.0            fastmap_1.2.0      
-#>  [7] janeaustenr_1.0.0   digest_0.6.38       rpart_4.1.24       
-#> [10] timechange_0.3.0    lifecycle_1.0.4     tokenizers_0.3.0   
-#> [13] survival_3.8-3      magrittr_2.0.4      compiler_4.5.1     
-#> [16] rlang_1.1.6         sass_0.4.10         tools_4.5.1        
-#> [19] utf8_1.2.6          yaml_2.3.10         tidytext_0.4.3     
-#> [22] data.table_1.17.8   knitr_1.50          askpass_1.2.1      
-#> [25] labeling_0.4.3      reticulate_1.44.0   RColorBrewer_1.1-3 
-#> [28] withr_3.0.2         purrr_1.2.0         nnet_7.3-20        
-#> [31] grid_4.5.1          embed_1.2.0         future_1.67.0      
-#> [34] ggplot2_4.0.0       globals_0.18.0      scales_1.4.0       
-#> [37] MASS_7.3-65         ggridges_0.5.7      cli_3.6.5          
-#> [40] rmarkdown_2.30      generics_0.1.4      umap_0.2.10.0      
-#> [43] RSpectra_0.16-2     future.apply_1.20.0 cachem_1.1.0       
-#> [46] splines_4.5.1       parallel_4.5.1      ggplotify_0.1.3    
-#> [49] BiocManager_1.30.25 yulab.utils_0.2.1   vctrs_0.6.5        
-#> [52] hardhat_1.4.2       Matrix_1.7-4        jsonlite_2.0.0     
-#> [55] gridGraphics_0.5-1  listenv_0.10.0      gower_1.0.2        
-#> [58] tidyr_1.3.1         jquerylib_0.1.4     recipes_1.3.1      
-#> [61] glue_1.8.0          parallelly_1.45.1   codetools_0.2-20   
-#> [64] uwot_0.2.4          RcppAnnoy_0.0.22    lubridate_1.9.4    
-#> [67] stringi_1.8.7       gtable_0.3.6        tibble_3.3.0       
-#> [70] pillar_1.11.1       rappdirs_0.3.3      htmltools_0.5.8.1  
-#> [73] openssl_2.3.4       ipred_0.9-15        lava_1.8.2         
-#> [76] R6_2.6.1            evaluate_1.0.5      lattice_0.22-6     
-#> [79] png_0.1-8           backports_1.5.0     SnowballC_0.7.1    
-#> [82] pheatmap_1.0.13     broom_1.0.10        renv_1.1.4         
-#> [85] bslib_0.9.0         class_7.3-23        Rcpp_1.1.0         
-#> [88] gridExtra_2.3       prodlim_2025.04.28  xfun_0.54          
-#> [91] fs_1.6.6            tidyheatmaps_0.2.1  pkgconfig_2.0.3
-
- - - - - - - - - - - diff --git a/inst/doc/differential_expression.html b/inst/doc/differential_expression.html deleted file mode 100644 index bbd2482..0000000 --- a/inst/doc/differential_expression.html +++ /dev/null @@ -1,568 +0,0 @@ - - - - - - - - - - - - - - -Differential Expression Analysis - - - - - - - - - - - - - - - - - - - - - - - - - - -

Differential Expression Analysis

- - - -

This vignette will guide you through the differential expression -analysis of your data. We will load HDAnalyzeR, load the example data -and metadata that come with the package and initialize the HDAnalyzeR -object.

-
-

Loading the Data

-
library(HDAnalyzeR)
-
-hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-
-
-

Running Differential Expression Analysis with limma

-

We will start by running a simple differential expression analysis -using the hd_de_limma() function. In this function we have -to state the variable of interest, the group of this variable that will -be the case, as well as the control(s). We will also correct for both -Sex and Age variables. After the analysis is -done, we will use hd_plot_volcano() to visualize the -results.

-

In the first example, we will run a differential expression analysis -for the AML case against the CLL control.

-
de_results <- hd_de_limma(hd_obj, 
-                          variable = "Disease", 
-                          case = "AML", 
-                          control = "CLL", 
-                          correct = c("Age", "Sex")) |> 
-  hd_plot_volcano()
-
-head(de_results$de_res)
-#> # A tibble: 6 × 10
-#>   Feature logFC   CI.L   CI.R AveExpr     t      P.Value adj.P.Val     B Disease
-#>   <chr>   <dbl>  <dbl>  <dbl>   <dbl> <dbl>        <dbl>     <dbl> <dbl> <chr>  
-#> 1 ADA      1.42  0.955  1.89    1.56   6.04 0.0000000250   2.50e-6  8.74 AML    
-#> 2 ADAM8   -1.23 -1.67  -0.795   1.74  -5.59 0.000000193    9.64e-6  6.77 AML    
-#> 3 AZU1     1.92  1.20   2.64    0.777  5.30 0.000000668    2.23e-5  5.57 AML    
-#> 4 ARID4B  -1.38 -1.91  -0.847   1.85  -5.15 0.00000132     3.30e-5  4.92 AML    
-#> 5 ARTN     1.08  0.597  1.55    0.804  4.46 0.0000227      3.84e-4  2.22 AML    
-#> 6 ANGPT1  -1.71 -2.47  -0.948   0.992 -4.45 0.0000230      3.84e-4  2.21 AML
-de_results$volcano_plot
-

-

We are able to state more control groups if we want to. We can also -change the correction for the variables as well as both the p-value and -logFC significance thresholds.

-
de_results <- hd_de_limma(hd_obj, 
-                              case = "AML", 
-                              control = c("CLL", "MYEL", "GLIOM"), 
-                              correct = "BMI") |> 
-  hd_plot_volcano(pval_lim = 0.01, logfc_lim  = 1)
-
-de_results$volcano_plot
-

-

If we do not set a control group, the function will compare the case -group against all other groups.

-
de_results <- hd_de_limma(hd_obj, case = "AML", correct = c("Age", "Sex")) |> 
-  hd_plot_volcano()
-
-de_results$volcano_plot
-

-
-
-

Customizing the Volcano Plot

-

We can customize the volcano plot further by adding a title and not -displaying the number of significant proteins. We can also change the -number of significant proteins that will be displayed with their names -in the plot.

-
de_results <- hd_de_limma(hd_obj, case = "AML", correct = c("Age", "Sex")) |> 
-  hd_plot_volcano(report_nproteins = FALSE, 
-                  title = "AML vs all other groups",
-                  top_up_prot = 3,
-                  top_down_prot = 1)
-
-de_results$volcano_plot
-

-
-
-

Running Differential Expression Analysis with t-test

-

Let’s move to another method. We will use the -hd_de_ttest() that performs a t-test for each variable. -This function takes similar inputs with hd_de_limma() but -it cannot correct for other variables like Sex and -Age.

-
de_results <- hd_de_ttest(hd_obj, case = "AML") |> 
-  hd_plot_volcano()
-
-de_results$volcano_plot
-

-
-
-

The case of Sex-specific Diseases

-

If we have diseases that are sex specific like Breast Cancer for -example, we should consider run the analysis only with samples of that -sex. We can easily integrate that into our pipeline using the -hd_filter() function. In that case, we would not be able to -correct for sex, as there will be only one sex “F” (female).

-
de_results <- hd_obj |> 
-  hd_filter(variable = "Sex", values = "F", flag = "k") |> 
-  hd_de_limma(case = "BRC", control = "AML", correct = "Age") |> 
-  hd_plot_volcano()
-#> VariableSexis categorical
-#> Filtering complete. Rows remaining:366
-
-de_results$volcano_plot
-

-
-
-

Running DE against other Variables

-
-

Other Categorical Variables

-

We could also run differential expression against another categorical -variable like Sex by changing the variable -argument.

-
de_results <- hd_de_limma(hd_obj, variable = "Sex", case = "F", correct = "Age") |> 
-  hd_plot_volcano(report_nproteins = FALSE, title = "Sex Comparison")
-
-de_results$volcano_plot
-

-
-
-

Continuous Variables

-

Moreover, we can also perform Differential Expression Analysis -against a continuous variable such as Age. This can be done -only with hd_de_limma()! We can also correct for -categorical and other continuous variables. In this case, no -case or control groups are needed.

-
de_results <- hd_de_limma(hd_obj, variable = "Age", case = NULL, correct = c("Sex", "BMI")) |> 
-  hd_plot_volcano(report_nproteins = FALSE, title = "DE against Age")
-
-de_results$volcano_plot
-

-
-
-
-

Summarizing the Results from Multiple Analysis

-

As a last step, we can summarize the results via -hd_plot_de_summary(). Let’s first run a differential -expression analysis for 4 different cases (1 vs 3).

-
res_aml <- hd_de_limma(hd_obj, case = "AML", control = c("CLL", "MYEL", "GLIOM"))
-
-res_cll <- hd_de_limma(hd_obj, case = "CLL", control = c("AML", "MYEL", "GLIOM"))
-
-res_myel <- hd_de_limma(hd_obj, case = "MYEL" , control = c("AML", "CLL", "GLIOM"))
-
-res_gliom <- hd_de_limma(hd_obj, case = "GLIOM" , control = c("AML", "CLL", "MYEL"))
-
-de_summary_res <- hd_plot_de_summary(list("AML" = res_aml, 
-                                          "CLL" = res_cll, 
-                                          "MYEL" = res_myel, 
-                                          "GLIOM" = res_gliom),
-                                     class_palette = "cancers12")
-
de_summary_res$de_barplot
-

-
de_summary_res$upset_plot_up
-

-
de_summary_res$upset_plot_down
-

-
-

📓 Remember that these data are a dummy-dataset with artificial data -and the results in this guide should not be interpreted as real results. -The purpose of this vignette is to show you how to use the package and -its functions.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats     graphics  grDevices datasets  utils     methods   base     
-#> 
-#> other attached packages:
-#> [1] HDAnalyzeR_1.0.0
-#> 
-#> loaded via a namespace (and not attached):
-#>  [1] gtable_0.3.6        limma_3.64.3        jsonlite_2.0.0     
-#>  [4] dplyr_1.1.4         compiler_4.5.1      BiocManager_1.30.25
-#>  [7] renv_1.1.4          Rcpp_1.1.0          tidyselect_1.2.1   
-#> [10] stringr_1.6.0       gridExtra_2.3       tidyr_1.3.1        
-#> [13] jquerylib_0.1.4     scales_1.4.0        yaml_2.3.10        
-#> [16] fastmap_1.2.0       statmod_1.5.1       plyr_1.8.9         
-#> [19] UpSetR_1.4.0        ggplot2_4.0.0       R6_2.6.1           
-#> [22] labeling_0.4.3      generics_0.1.4      knitr_1.50         
-#> [25] MASS_7.3-65         ggrepel_0.9.6       tibble_3.3.0       
-#> [28] RColorBrewer_1.1-3  bslib_0.9.0         pillar_1.11.1      
-#> [31] rlang_1.1.6         utf8_1.2.6          stringi_1.8.7      
-#> [34] cachem_1.1.0        xfun_0.54           S7_0.2.0           
-#> [37] sass_0.4.10         cli_3.6.5           withr_3.0.2        
-#> [40] magrittr_2.0.4      digest_0.6.38       grid_4.5.1         
-#> [43] lifecycle_1.0.4     vctrs_0.6.5         evaluate_1.0.5     
-#> [46] glue_1.8.0          farver_2.1.2        rmarkdown_2.30     
-#> [49] purrr_1.2.0         tools_4.5.1         pkgconfig_2.0.3    
-#> [52] htmltools_0.5.8.1
-
- - - - - - - - - - - diff --git a/inst/doc/post_analysis.html b/inst/doc/post_analysis.html deleted file mode 100644 index 8c583f1..0000000 --- a/inst/doc/post_analysis.html +++ /dev/null @@ -1,540 +0,0 @@ - - - - - - - - - - - - - - -Post Analysis: Pathway Enrichment - - - - - - - - - - - - - - - - - - - - - - - - - - -

Post Analysis: Pathway Enrichment

- - - -

This vignette will guide you through the post analysis of the results -obtained from the HDAnalyzeR pipeline. The pathway enrichment analysis -is performed using the Gene Ontology, KEGG and Reactome databases from -clusterProfiler and ReactomePA packages -respectively.

-

If you want to learn more about ORA and GSEA, please refer to the -following publications:

- -
-

📓 Remember that these data are a dummy-dataset with artificial data -and the results in this guide should not be interpreted as real results. -This is why we are using extremely large p-value cutoffs in this case -that should not be used in real data.

-
-
-

Loading the Data

-

We will load HDAnalyzeR and dplyr, load the example data and metadata -that come with the package and initialize the HDAnalyzeR object.

-
library(HDAnalyzeR)
-library(dplyr)
-
-hd_obj <- hd_initialize(dat = example_data, 
-                        metadata = example_metadata, 
-                        is_wide = FALSE, 
-                        sample_id = "DAid",
-                        var_name = "Assay",
-                        value_name = "NPX")
-

For the Over Representation Analysis we are going to use a list of -differentially expressed proteins. In this example we are going to use -the up-regulated proteins. We could also use the features list from the -classification models or even run both and get the intersect as it is -done in the Get Started guide.

-
de_res <- hd_de_limma(hd_obj, case = "AML")
-
-
-

Over Representation Analysis

-

First, we will perform an Over Representation Analysis (ORA) using -the Gene Ontology database and the BP ontology. We will use the -hd_ora() and hd_plot_ora() functions to run -the analysis and plot the results respectively.

-
proteins <- de_res$de_res |> 
-  filter(logFC > 0 & adj.P.Val < 0.05) |> 
-  pull(Feature)
-
-enrichment <- hd_ora(proteins, database = "GO", ontology = "BP")
-
-enrichment_plots <- hd_plot_ora(enrichment)
-
-enrichment_plots$dotplot
-

-
enrichment_plots$treeplot
-

-
enrichment_plots$cnetplot
-

-

Let’s change the database and the p-value threshold.

-
enrichment  <- hd_ora(proteins, database = "Reactome", pval_lim = 0.2)
-
-enrichment_plots <- hd_plot_ora(enrichment)
-
-enrichment_plots$dotplot
-

-
enrichment_plots$treeplot
-

-
enrichment_plots$cnetplot
-

-
-
-

Gene Set Enrichment Analysis

-

We can also run a Gene Set Enrichment Analysis (GSEA) using the -hd_gsea() and hd_plot_gsea functions. The -hd_plot_gsea() function will plot the results.

-
-

⚠️ In this case, the function requires strictly differential -expression results, so a ranked list of proteins is derived based on the -ranked_by argument.

-
-
enrichment <- hd_gsea(de_res, database = "GO", ontology = "BP", pval_lim = 0.55)
-
-enrichment_plots <- hd_plot_gsea(enrichment)
-
-enrichment_plots$dotplot
-

-
enrichment_plots$gseaplot
-

-
enrichment_plots$cnetplot
-

-
enrichment_plots$ridgeplot
-

-

We can also change the ranking variable to the product of logFC and --log(adjusted p value) instead of the default logFC by changing the -ranked_by argument to “both”. We could also use other -variables such as p-value or any other variable in the DE results. -However, you should use as ranking a variable that has some form of -biological relevance of the variable.

-
enrichment <- hd_gsea(de_res, 
-                      database = "GO", 
-                      ontology = "BP", 
-                      pval_lim = 0.9, 
-                      ranked_by = "both")
-
-enrichment_plots <- hd_plot_gsea(enrichment)
-enrichment_plots$cnetplot
-

-
-

📓 Remember once again that these data are a dummy-dataset with -artificial data and the results in this guide should not be interpreted -as real results. The purpose of this vignette is to show you how to use -the package and its functions.

-
-
sessionInfo()
-#> R version 4.5.1 (2025-06-13)
-#> Platform: aarch64-apple-darwin20
-#> Running under: macOS Tahoe 26.1
-#> 
-#> Matrix products: default
-#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
-#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
-#> 
-#> locale:
-#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
-#> 
-#> time zone: Europe/Stockholm
-#> tzcode source: internal
-#> 
-#> attached base packages:
-#> [1] stats4    stats     graphics  grDevices datasets  utils     methods  
-#> [8] base     
-#> 
-#> other attached packages:
-#> [1] org.Hs.eg.db_3.21.0  AnnotationDbi_1.70.0 IRanges_2.42.0      
-#> [4] S4Vectors_0.46.0     Biobase_2.68.0       BiocGenerics_0.54.1 
-#> [7] generics_0.1.4       dplyr_1.1.4          HDAnalyzeR_1.0.0    
-#> 
-#> loaded via a namespace (and not attached):
-#>   [1] DBI_1.2.3               gson_0.1.0              gridExtra_2.3          
-#>   [4] rlang_1.1.6             magrittr_2.0.4          DOSE_4.2.0             
-#>   [7] ggridges_0.5.7          compiler_4.5.1          RSQLite_2.4.4          
-#>  [10] reactome.db_1.92.0      png_0.1-8               vctrs_0.6.5            
-#>  [13] reshape2_1.4.5          stringr_1.6.0           pkgconfig_2.0.3        
-#>  [16] crayon_1.5.3            fastmap_1.2.0           XVector_0.48.0         
-#>  [19] ggraph_2.2.2            labeling_0.4.3          rmarkdown_2.30         
-#>  [22] enrichplot_1.28.4       graph_1.86.0            UCSC.utils_1.4.0       
-#>  [25] purrr_1.2.0             bit_4.6.0               xfun_0.54              
-#>  [28] cachem_1.1.0            graphite_1.54.0         aplot_0.2.9            
-#>  [31] GenomeInfoDb_1.44.3     jsonlite_2.0.0          blob_1.2.4             
-#>  [34] tweenr_2.0.3            BiocParallel_1.42.2     parallel_4.5.1         
-#>  [37] R6_2.6.1                bslib_0.9.0             stringi_1.8.7          
-#>  [40] RColorBrewer_1.1-3      limma_3.64.3            jquerylib_0.1.4        
-#>  [43] GOSemSim_2.34.0         Rcpp_1.1.0              knitr_1.50             
-#>  [46] ggtangle_0.0.8          R.utils_2.13.0          Matrix_1.7-4           
-#>  [49] splines_4.5.1           igraph_2.2.1            tidyselect_1.2.1       
-#>  [52] viridis_0.6.5           qvalue_2.40.0           yaml_2.3.10            
-#>  [55] codetools_0.2-20        lattice_0.22-6          tibble_3.3.0           
-#>  [58] plyr_1.8.9              treeio_1.32.0           withr_3.0.2            
-#>  [61] KEGGREST_1.48.1         S7_0.2.0                evaluate_1.0.5         
-#>  [64] gridGraphics_0.5-1      polyclip_1.10-7         Biostrings_2.76.0      
-#>  [67] pillar_1.11.1           BiocManager_1.30.25     ggtree_3.16.3          
-#>  [70] renv_1.1.4              clusterProfiler_4.16.0  ggfun_0.2.0            
-#>  [73] ggplot2_4.0.0           scales_1.4.0            tidytree_0.4.6         
-#>  [76] glue_1.8.0              lazyeval_0.2.2          tools_4.5.1            
-#>  [79] ggnewscale_0.5.2        data.table_1.17.8       fgsea_1.34.2           
-#>  [82] graphlayouts_1.2.2      fs_1.6.6                tidygraph_1.3.1        
-#>  [85] fastmatch_1.1-6         cowplot_1.2.0           grid_4.5.1             
-#>  [88] tidyr_1.3.1             ape_5.8-1               colorspace_2.1-2       
-#>  [91] nlme_3.1-168            GenomeInfoDbData_1.2.14 patchwork_1.3.2        
-#>  [94] ggforce_0.5.0           cli_3.6.5               rappdirs_0.3.3         
-#>  [97] viridisLite_0.4.2       ReactomePA_1.52.0       gtable_0.3.6           
-#> [100] R.methodsS3_1.8.2       yulab.utils_0.2.1       sass_0.4.10            
-#> [103] digest_0.6.38           ggrepel_0.9.6           ggplotify_0.1.3        
-#> [106] farver_2.1.2            memoise_2.0.1           htmltools_0.5.8.1      
-#> [109] R.oo_1.27.1             lifecycle_1.0.4         httr_1.4.7             
-#> [112] GO.db_3.21.0            statmod_1.5.1           bit64_4.6.0-1          
-#> [115] MASS_7.3-65
-
- - - - - - - - - - - diff --git a/man/hd_plot_feature_heatmap.Rd b/man/hd_plot_feature_heatmap.Rd index 995ee68..1f1684f 100644 --- a/man/hd_plot_feature_heatmap.Rd +++ b/man/hd_plot_feature_heatmap.Rd @@ -30,7 +30,7 @@ The summary heatmap of the combined differential expression and classification m \code{hd_plot_feature_heatmap} plots a summary bubble-heatmap of the combined differential expression and classification model results. The heatmap shows the log2 fold change and adjusted p-value of the differential expression results, and the scaled importance -and sign of the classification model results. The heatmap is ordered and the selected +of the classification model results. The heatmap is ordered and the selected assays are based on the specified control group. * are added to the features that pass the p-value threshold. } diff --git a/man/variable_imp.Rd b/man/variable_imp.Rd index debb1c0..2792725 100644 --- a/man/variable_imp.Rd +++ b/man/variable_imp.Rd @@ -13,6 +13,7 @@ variable_imp( y_labels = FALSE, title = c("accuracy", "sensitivity", "specificity", "auc", "features", "top-features"), verbose = TRUE, + engine = "glmnet", seed = 123 ) } @@ -29,6 +30,8 @@ variable_imp( \item{title}{Vector of title elements to include in the plot.} +\item{engine}{Either glmnet, rf or lr. If glmnet and the model type is multi-class then the variable importance is calculated per class.} + \item{seed}{Seed for reproducibility. Default is 123.} } \value{ From eb5a5f45f9bffcca66814c2bd09303bceb7baa3b Mon Sep 17 00:00:00 2001 From: Konstantinos Antonopoulos Date: Tue, 9 Dec 2025 16:31:26 +0100 Subject: [PATCH 3/4] Brought back sign in glmnet classification and also fix it for rf models by making stating that all signs there are positive which is correct conceptually for this model type --- R/classification_models.R | 8 ++++++-- R/visualize.R | 8 +++++--- man/hd_plot_feature_heatmap.Rd | 2 +- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/R/classification_models.R b/R/classification_models.R index 18073f2..1c6df3a 100644 --- a/R/classification_models.R +++ b/R/classification_models.R @@ -1077,6 +1077,7 @@ variable_imp <- function(dat, dplyr::filter(!!rlang::sym("term") != "(Intercept)") |> dplyr::select(-dplyr::any_of(c("penalty"))) |> dplyr::mutate( + Sign = ifelse(!!rlang::sym("estimate") > 0, "POS", "NEG"), Feature = !!rlang::sym("term"), Importance = abs(!!rlang::sym("estimate")) ) |> @@ -1093,8 +1094,11 @@ variable_imp <- function(dat, dplyr::ungroup() } else { features <- features |> - dplyr::mutate(Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance"))) |> - dplyr::mutate(Feature = forcats::fct_reorder(!!rlang::sym("Feature"), !!rlang::sym("Importance"))) + dplyr::mutate( + Sign = "POS", + Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance")), + Feature = forcats::fct_reorder(!!rlang::sym("Feature"), !!rlang::sym("Importance")) + ) } } else { diff --git a/R/visualize.R b/R/visualize.R index e6446f7..c090456 100644 --- a/R/visualize.R +++ b/R/visualize.R @@ -320,7 +320,7 @@ hd_plot_regression <- function(dat, #' `hd_plot_feature_heatmap` plots a summary bubble-heatmap of the combined differential #' expression and classification model results. The heatmap shows the log2 fold change #' and adjusted p-value of the differential expression results, and the scaled importance -#' of the classification model results. The heatmap is ordered and the selected +#' and sign of the classification model results. The heatmap is ordered and the selected #' assays are based on the specified control group. * are added to the features that pass #' the p-value threshold. #' @@ -414,7 +414,8 @@ hd_plot_feature_heatmap <- function(de_results, res_model <- model_results[[i]][["features"]] |> dplyr::filter(!!rlang::sym("Feature") %in% assays) |> dplyr::select(!!rlang::sym("Feature"), - !!rlang::sym("Scaled_Importance")) + !!rlang::sym("Scaled_Importance"), + !!rlang::sym("Sign")) res_combined <- res_de |> dplyr::left_join(res_model, by = c("Feature")) |> @@ -433,10 +434,11 @@ hd_plot_feature_heatmap <- function(de_results, ) |> ggplot2::ggplot(ggplot2::aes(x = !!rlang::sym("Feature"), y = !!rlang::sym("control_group"))) + ggplot2::geom_tile(ggplot2::aes(fill = !!rlang::sym("logFC")), color = "white") + - ggplot2::geom_point(ggplot2::aes(size = !!rlang::sym("Scaled_Importance"))) + + ggplot2::geom_point(ggplot2::aes(size = !!rlang::sym("Scaled_Importance"), color = !!rlang::sym("Sign"))) + ggplot2::geom_point(ggplot2::aes(size = !!rlang::sym("Scaled_Importance")), shape = 1, colour = "black") + ggplot2::geom_text(ggplot2::aes(label = ifelse(!!rlang::sym("adj.P.Val") < pval_lim, "*", "")), color = "black", size = 3) + ggplot2::scale_fill_gradient2(low = "#317EC2", mid = "white", high = "#C03830", midpoint = 0, name = "Log2 FC") + + ggplot2::scale_color_manual(values = c("NEG" = "#317EC2", "POS" = "#C03830"), name = "Sign", na.translate = FALSE) + ggplot2::scale_size(name = "Importance") + ggplot2::labs(x = "Feature", y = "Control Group") + ggplot2::theme_minimal() + diff --git a/man/hd_plot_feature_heatmap.Rd b/man/hd_plot_feature_heatmap.Rd index 1f1684f..995ee68 100644 --- a/man/hd_plot_feature_heatmap.Rd +++ b/man/hd_plot_feature_heatmap.Rd @@ -30,7 +30,7 @@ The summary heatmap of the combined differential expression and classification m \code{hd_plot_feature_heatmap} plots a summary bubble-heatmap of the combined differential expression and classification model results. The heatmap shows the log2 fold change and adjusted p-value of the differential expression results, and the scaled importance -of the classification model results. The heatmap is ordered and the selected +and sign of the classification model results. The heatmap is ordered and the selected assays are based on the specified control group. * are added to the features that pass the p-value threshold. } From 1b8b31d71d060b936ddfd7da295e0cdc7ea814f3 Mon Sep 17 00:00:00 2001 From: Konstantinos Antonopoulos Date: Tue, 9 Dec 2025 22:36:34 +0100 Subject: [PATCH 4/4] Update to handle and return features correctly in all cases --- R/classification_models.R | 37 +++++++++++++++++++++++++++---------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/R/classification_models.R b/R/classification_models.R index 1c6df3a..7cbafd0 100644 --- a/R/classification_models.R +++ b/R/classification_models.R @@ -1086,21 +1086,22 @@ variable_imp <- function(dat, # Scale importance if (model_type == "multi_class") { features <- features |> - dplyr::group_by(class) |> + dplyr::group_by(!!rlang::sym("class")) |> dplyr::mutate( Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance")), - Feature_plot = tidytext::reorder_within(!!rlang::sym("Feature"), !!rlang::sym("Importance"), !!rlang::sym("class")) + Feature_plot = tidytext::reorder_within(!!rlang::sym("Feature"), !!rlang::sym("Importance"), !!rlang::sym("class")), + Class = !!rlang::sym("class") ) |> - dplyr::ungroup() + dplyr::ungroup() |> + dplyr::select(dplyr::any_of(c("Class", "Feature", "Importance", "Sign", "Scaled_Importance", "Feature_plot"))) } else { features <- features |> dplyr::mutate( - Sign = "POS", Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance")), Feature = forcats::fct_reorder(!!rlang::sym("Feature"), !!rlang::sym("Importance")) - ) + ) |> + dplyr::select(dplyr::any_of(c("Feature", "Importance", "Sign", "Scaled_Importance"))) } - } else { features <- final |> workflows::extract_fit_parsnip() |> @@ -1112,8 +1113,16 @@ variable_imp <- function(dat, Variable = forcats::fct_reorder(!!rlang::sym("Variable"), !!rlang::sym("Importance"))) |> dplyr::arrange(dplyr::desc(!!rlang::sym("Importance"))) |> # Min max scaling with min = 0 always and max = 1 - dplyr::mutate(Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance"))) |> + dplyr::mutate( + Scaled_Importance = !!rlang::sym("Importance") / max(!!rlang::sym("Importance")) + ) |> dplyr::rename(Feature = !!rlang::sym("Variable")) + + if (engine == "rf") { + features <- features |> + dplyr::mutate(Sign = "POS") |> + dplyr::relocate(!!rlang::sym("Sign"), .after = !!rlang::sym("Importance")) + } } if (model_type == "binary_class") { @@ -1153,7 +1162,7 @@ variable_imp <- function(dat, } else if (!is.null(palette)) { pal <- palette } else { - classes <- unique(features$class) + classes <- unique(features[["Class"]]) pal <- rep("#883268", length(classes)) } @@ -1184,8 +1193,8 @@ variable_imp <- function(dat, var_imp_plot <- features |> dplyr::filter(!!rlang::sym("Scaled_Importance") > 0) |> ggplot2::ggplot(ggplot2::aes(x = !!rlang::sym("Feature_plot"), y = !!rlang::sym("Scaled_Importance"))) + - ggplot2::geom_col(ggplot2::aes(fill = ifelse(!!rlang::sym("Scaled_Importance") > 0.5, class, NA))) + - ggplot2::facet_wrap(~ class, scales = "free_y") + + ggplot2::geom_col(ggplot2::aes(fill = ifelse(!!rlang::sym("Scaled_Importance") > 0.5, !!rlang::sym("Class"), NA))) + + ggplot2::facet_wrap(~ Class, scales = "free_y") + tidytext::scale_x_reordered() + # cleans axis labels ggplot2::coord_flip() + # horizontal bars ggplot2::scale_fill_manual(values = pal, na.value = "grey80") + @@ -1196,6 +1205,14 @@ variable_imp <- function(dat, axis.text.x = ggplot2::element_text(hjust = 0.5), panel.spacing = ggplot2::unit(1, "lines") ) + + features <- features |> + dplyr::select(-dplyr::any_of(c("Feature_plot"))) |> + dplyr::arrange( + !!rlang::sym("Class"), + dplyr::desc(!!rlang::sym("Scaled_Importance")) + ) + } else { var_imp_plot <- features |> dplyr::filter(!!rlang::sym("Scaled_Importance") > 0) |>