diff --git a/.Rbuildignore b/.Rbuildignore
index 5d97a0bf..d713df39 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -20,4 +20,5 @@
 ^\\.git$
 ^codecov\.yml$
 ^vignettes/articles$
-^data_aux$
\ No newline at end of file
+^data_aux$ 
+^python(/.*)?$
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index ea55e2da..32e172ac 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,4 @@ markeR.Rproj
 inst/doc
 data_aux
 markeR.Rcheck 
+/python/.venv
diff --git a/DESCRIPTION b/DESCRIPTION
index 1d73225d..17cce5fd 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: markeR
 Title: An R Toolkit for Evaluating Gene Signatures as Phenotypic Markers
-Version: 1.1.0
+Version: 1.1.2
 Authors@R: 
     c(
       person("Rita", "Martins-Silva", 
@@ -25,7 +25,7 @@ Encoding: UTF-8
 Language: en-GB
 LazyData: false
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Additional_repositories: https://bioconductor.org/packages/release/bioc
 Imports:
     circlize,
diff --git a/NEWS.md b/NEWS.md
index 538bb1eb..4d66e94e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,28 @@
+# markeR 1.1.2 (12 Mar, 2026) 
+
+## Minor Changes
+- Moved Python bridge scripts from `inst/python/` to a top-level `python/` 
+  directory, as these are supplementary scripts not part of the R package itself.
+- Added `requirements.txt` to the `python/` directory listing all needed 
+  Python dependencies (`rpy2`, `pandas`, `numpy`, and optionally 
+  `ipython` and `jupyter`) for easier environment setup. 
+- Removed redundant code snippets from the Python bridge scripts.
+
+# markeR 1.1.1 (11 Mar, 2026) 
+ 
+  - Added `p.adjust.method` parameter across all functions performing or
+  depending on multiple testing correction, allowing users to specify
+  any correction method supported by `stats::p.adjust()`, beyond the default 
+  Benjamini-Hochberg FDR.
+- Added Python bridge scripts in `inst/python/` for users who wish to call
+  markeR from a Python environment via `rpy2`. Includes a tutorial workflow
+  script and a generic command-line wrapper capable of invoking any exported
+  markeR function. See `inst/python/README.md` for installation and usage.
+
+# markeR 1.0.0 (31 Oct, 2025)
+
+- Official Bioconductor Release.
+
 # markeR 0.99.5 (17 Sep, 2025)
 
 - Minor fix in `.onAttach()` to avoid errors when checking `ggplot2` version and ensure the startup warning works correctly.
diff --git a/R/FPR_Simulation.R b/R/FPR_Simulation.R
index 02c760b0..126e44ad 100644
--- a/R/FPR_Simulation.R
+++ b/R/FPR_Simulation.R
@@ -44,7 +44,11 @@ utils::globalVariables(c( "cohen", "method", "contrast" ))
 #'   grid layout. If `NULL`, layout is auto-calculated.
 #' @param nrow Integer. Number of rows for arranging signature plots in a grid
 #'   layout. If `NULL`, layout is auto-calculated.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
 #'
 #' @return Invisibly returns a list containing:
 #'   \describe{
@@ -111,7 +115,7 @@ FPR_Simulation <- function(data, metadata, original_signatures, Variable,
                            gene_list = NULL, number_of_sims=100, title=NULL,
                            widthTitle = 30, titlesize = 12,  pointSize = 2,
                            labsize = 10,mode = c( "none","simple","medium","extensive"),
-                           ColorValues=NULL, ncol=NULL, nrow=NULL) {
+                           ColorValues=NULL, ncol=NULL, nrow=NULL, p.adjust.method="BH") {
   data <- as.data.frame(data) # Ensure data is a data frame
   if (is.null(gene_list)) gene_list <- row.names(data)
 
@@ -128,7 +132,7 @@ FPR_Simulation <- function(data, metadata, original_signatures, Variable,
       results <- suppressMessages(CohenF_allConditions(data = data,
                                                        metadata = metadata,
                                                        gene_sets = original_signatures,
-                                                       variable = Variable ))
+                                                       variable = Variable, p.adjust.method = p.adjust.method ))
       cohentype <- "f"
 
     } else {
@@ -138,7 +142,7 @@ FPR_Simulation <- function(data, metadata, original_signatures, Variable,
         results <- suppressMessages(CohenF_allConditions(data = data,
                                                          metadata = metadata,
                                                          gene_sets = original_signatures,
-                                                         variable = Variable ))
+                                                         variable = Variable, p.adjust.method = p.adjust.method ))
         cohentype <- "f"
 
       } else {
@@ -146,7 +150,7 @@ FPR_Simulation <- function(data, metadata, original_signatures, Variable,
         results <- suppressMessages(CohenD_allConditions(data = data,
                                                          metadata = metadata,
                                                          gene_sets = original_signatures,
-                                                         variable = Variable, mode = mode))
+                                                         variable = Variable, mode = mode, p.adjust.method = p.adjust.method))
         cohentype <- "d"
 
       }
@@ -219,14 +223,14 @@ FPR_Simulation <- function(data, metadata, original_signatures, Variable,
         metadata = metadata,
         gene_sets = simulatedsigs,
         variable = Variable,
-        mode = mode
+        mode = mode, p.adjust.method = p.adjust.method
       ))
     } else {
       results2 <- suppressMessages(CohenF_allConditions(
         data = data,
         metadata = metadata,
         gene_sets = simulatedsigs,
-        variable = Variable
+        variable = Variable, p.adjust.method = p.adjust.method
       ))
     }
 
diff --git a/R/GSEA_VariableAssociation.R b/R/GSEA_VariableAssociation.R
index 1547a9ed..159d8156 100644
--- a/R/GSEA_VariableAssociation.R
+++ b/R/GSEA_VariableAssociation.R
@@ -57,7 +57,12 @@
 #'   removed before analysis, leading to a loss of data to be fitted in the
 #'   model.
 #' @param printplt Boolean specifying if plot is to be printed. Default: `TRUE`.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @return A list with two elements:
 #'   - `data`: A data frame containing the GSEA results, including normalized
 #'   enrichment scores (NES), adjusted p-values, and contrasts.
@@ -70,7 +75,8 @@ GSEA_VariableAssociation <- function(data, metadata, cols, stat=NULL,
                                      signif_color = "red", saturation_value=NULL,
                                      sig_threshold = 0.05, widthlabels=18,
                                      labsize=10, titlesize=14, pointSize=5,
-                                     ignore_NAs = FALSE, printplt =TRUE) {
+                                     ignore_NAs = FALSE, printplt =TRUE,
+                                     p.adjust.method = "BH") {
   data <- as.data.frame(data) # Ensure data is a data frame
   mode <- match.arg(mode)
   metadata <- metadata[, cols %in% colnames(metadata), drop = FALSE]
@@ -129,7 +135,7 @@ GSEA_VariableAssociation <- function(data, metadata, cols, stat=NULL,
   combined_results$Contrast <- cont_vec
 
   # correct adjusted p value to correct for multiple testing for the contrasts?
-  combined_results$padj <- stats::p.adjust(combined_results$padj, method = "BH")
+  combined_results$padj <- stats::p.adjust(combined_results$padj, method = p.adjust.method)
 
 
   combined_results_toreturn <- combined_results
diff --git a/R/Heatmap_Cohen.R b/R/Heatmap_Cohen.R
index 241ebd4c..89aac9d3 100644
--- a/R/Heatmap_Cohen.R
+++ b/R/Heatmap_Cohen.R
@@ -33,6 +33,7 @@
 #' @param ColorValues A character vector specifying the colors for the gradient
 #'   fill in the heatmaps. Default is \code{c("#F9F4AE", "#B44141")}.
 #' @param title Title for the grid of plots.
+#' 
 #' @return A list with two elements:
 #' \describe{
 #'   \item{plt}{A combined heatmap arranged in a grid using \code{ggpubr::ggarrange}.}
@@ -61,7 +62,7 @@
 #' @keywords internal
 Heatmap_Cohen <- function(cohenlist, nrow = NULL, ncol = NULL, limits = NULL,
                           widthTitle = 22, titlesize = 12, ColorValues = NULL,
-                          title=NULL ) {
+                          title=NULL) {
 
   cohentype <- ifelse("CohenD" %in% names(cohenlist[[1]]), "d",
                       ifelse("CohenF" %in% names(cohenlist[[1]]), "f", NULL))
@@ -201,7 +202,12 @@ Heatmap_Cohen <- function(cohenlist, nrow = NULL, ncol = NULL, limits = NULL,
 #' groups.
 #' - `"extensive"`: All possible groupwise contrasts, ensuring balance in the
 #' number of terms on each side.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @return A named list where each element corresponds to a gene signature. Each
 #'   signature element is a list with three components:
 #' \describe{
@@ -227,7 +233,7 @@ Heatmap_Cohen <- function(cohenlist, nrow = NULL, ncol = NULL, limits = NULL,
 #'
 #' @keywords internal
 CohenD_allConditions <- function(data, metadata, gene_sets, variable,
-                                 mode = c("simple","medium","extensive")) {
+                                 mode = c("simple","medium","extensive"), p.adjust.method = "BH") {
   data <- as.data.frame(data) # Ensure data is a data frame
   # Step 1: Check if variable exists in metadata
   if (!variable %in% colnames(metadata)) {
@@ -292,8 +298,8 @@ CohenD_allConditions <- function(data, metadata, gene_sets, variable,
     }
   }
 
-  # Step 2: Apply BH correction within each method
-  all_padj <- lapply(all_pvalues, function(pvals) stats::p.adjust(pvals, method = "BH"))
+  # Step 2: Apply  correction within each method
+  all_padj <- lapply(all_pvalues, function(pvals) stats::p.adjust(pvals, method = p.adjust.method))
 
   # Step 3: Store corrected p-values back into result_list
   index_tracker <- list()  # Track index position for each method
@@ -495,6 +501,12 @@ flatten_results <- function(nested_list) {
 #'   downregulated).
 #' @param variable A string specifying the categorical variable in
 #'   \code{metadata} used to model the gene signature scores.
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @return A named list where each element corresponds to a gene signature. Each
 #'   signature element is a list with three components:
 #' \describe{
@@ -508,7 +520,7 @@ flatten_results <- function(nested_list) {
 #' }
 #'
 #' @keywords internal
-CohenF_allConditions <- function(data, metadata, gene_sets, variable ) {
+CohenF_allConditions <- function(data, metadata, gene_sets, variable, p.adjust.method = "BH" ) {
   data <- as.data.frame(data) # Ensure data is a data frame
   # Step 1: Check if variable exists in metadata
   if (!variable %in% colnames(metadata)) {
@@ -579,8 +591,8 @@ CohenF_allConditions <- function(data, metadata, gene_sets, variable ) {
     }
   }
 
-  # Step 2: Apply BH correction within each method
-  all_padj <- lapply(all_pvalues, function(pvals) stats::p.adjust(pvals, method = "BH"))
+  # Step 2: Apply   correction within each method
+  all_padj <- lapply(all_pvalues, function(pvals) stats::p.adjust(pvals, method = p.adjust.method))
 
   # Step 3: Store corrected p-values back into result_list
   index_tracker <- list()  # Track index position for each method
diff --git a/R/PlotScores.R b/R/PlotScores.R
index 8ab0ad4f..0dc36cf2 100644
--- a/R/PlotScores.R
+++ b/R/PlotScores.R
@@ -114,7 +114,12 @@ utils::globalVariables(c("score"))
 #' @param cor Correlation method for numeric variables. One of `"pearson"`
 #'   (default), `"spearman"`, or `"kendall"`. Only applies when the variable is
 #'   numeric and `method != "all"`.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. Only if `method == "all"`.
+#'   
 #' @return Depending on `method`:
 #'
 #'   If `method = "all"`, returns a list with `heatmap` and `volcano` ggplot objects.
@@ -241,7 +246,8 @@ PlotScores <- function(data, metadata, gene_sets,
                        cond_cohend = NULL, pvalcalc = FALSE,
                        mode = c("simple","medium","extensive"),
                        widthlegend=22, sig_threshold=0.05, cohen_threshold=0.5,
-                       colorPalette="Set3", cor=c("pearson","spearman","kendall")) {
+                       colorPalette="Set3", cor=c("pearson","spearman","kendall"),
+                       p.adjust.method="BH") {
   data <- as.data.frame(data) # Ensure data is a data frame
   method <- match.arg(method)
   mode <- match.arg(mode)
@@ -254,13 +260,13 @@ PlotScores <- function(data, metadata, gene_sets,
     if (type =="Numeric"){
 
       cohenlist <- CohenF_allConditions(data = data, metadata = metadata,
-                                        gene_sets = gene_sets, variable = Variable )
+                                        gene_sets = gene_sets, variable = Variable, p.adjust.method = p.adjust.method )
 
     } else {
 
       cohenlist <- CohenD_allConditions(data = data, metadata = metadata,
                                         gene_sets = gene_sets, variable = Variable,
-                                        mode = mode)
+                                        mode = mode, p.adjust.method = p.adjust.method )
 
     }
 
diff --git a/R/Score_VariableAssociation.R b/R/Score_VariableAssociation.R
index 7c52e0cc..8c38b616 100644
--- a/R/Score_VariableAssociation.R
+++ b/R/Score_VariableAssociation.R
@@ -161,7 +161,12 @@ create_contrast_column <- function(metadata, variable_name, contrast) {
 #' @param color_palette A string specifying the color palette for discrete
 #' variables. Default: `"Set2"`.
 #' @param printplt Boolean specifying if plot is to be printed. Default: `TRUE`.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @return A list with:
 #'   - `Overall`: Data frame of effect sizes and p-values for each contrasted
 #'   phenotypic variable.
@@ -207,7 +212,8 @@ Score_VariableAssociation <- function(data,
                                       discrete_colors=NULL,
                                       continuous_color = "#8C6D03",
                                       color_palette = "Set2",
-                                      printplt =TRUE){
+                                      printplt =TRUE, 
+                                      p.adjust.method = "BH"){
   method <- match.arg(method)  # Validate method input
   mode <- match.arg(mode)
   data <- as.data.frame(data) # Ensure data is a data frame
@@ -296,8 +302,7 @@ Score_VariableAssociation <- function(data,
   # Would happen if we have only numeric variables
   if (nrow(df_results_contrast)!=0){
 
-    df_results_contrast$padj <- stats::p.adjust(df_results_contrast$PValue,
-                                         method = "BH")
+    df_results_contrast$padj <- stats::p.adjust(df_results_contrast$PValue,  method = p.adjust.method)
 
     if(is.null(saturation_value)){
       if (min(df_results_contrast$padj)>sig_threshold){
diff --git a/R/VariableAssociation.R b/R/VariableAssociation.R
index 35c423d2..a321be0f 100644
--- a/R/VariableAssociation.R
+++ b/R/VariableAssociation.R
@@ -96,7 +96,12 @@ identify_variable_type <- function(df, cols = NULL) {
 #' @param categorical_multi The statistical test for multi-level categorical
 #' variables.
 #'   Options: `"anova"` (default) or `"kruskal-wallis"`.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @return A named list (one entry per variable being analysed) where each
 #' element is a data frame with:
 #'   - **Metric**: The test statistic (correlation coefficient, t-statistic,
@@ -130,7 +135,7 @@ identify_variable_type <- function(df, cols = NULL) {
 compute_stat_tests <- function(df, target_var, cols = NULL,
                                numeric = "pearson",
                                categorical_bin = "t.test",
-                               categorical_multi = "anova") {
+                               categorical_multi = "anova", p.adjust.method="BH") {
 
   # Ensure only one method is selected per variable type
   if (length(numeric) > 1 | length(categorical_bin) > 1 |
@@ -202,7 +207,7 @@ compute_stat_tests <- function(df, target_var, cols = NULL,
     # scientific notation
     test_df$metric <- formatC(test_df$metric, format = "e", digits = 2)
     # correct for multiple testing per variable
-    test_df$p_value <- stats::p.adjust(test_df$p_value, method = "BH")
+    test_df$p_value <- stats::p.adjust(test_df$p_value,  method = p.adjust.method)
     test_df$p_value <- formatC(test_df$p_value, format = "e", digits = 3)
 
 
@@ -267,7 +272,12 @@ compute_stat_tests <- function(df, target_var, cols = NULL,
 #' (`"B"` or `"t"`). Auto-detected if `NULL`.
 #' @param ignore_NAs (GSEA only) Logical. If `TRUE`, rows with NA metadata are
 #' removed. Default: `FALSE`.
-#'
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @return A list with method-specific results and ggplot2-based visualizations:
 #'
 #' **For score-based methods (`logmedian`, `ssGSEA`, `ranking`):**
@@ -355,7 +365,8 @@ VariableAssociation <- function(method = c("ssGSEA", "logmedian",
                                 discrete_colors = NULL,
                                 continuous_color = "#8C6D03",
                                 color_palette = "Set2",
-                                printplt = TRUE) {
+                                printplt = TRUE, 
+                                p.adjust.method = "BH") {
   method <- match.arg(method)
   mode <- match.arg(mode)
   data <- as.data.frame(data) # Ensure data is a data frame
@@ -375,7 +386,8 @@ VariableAssociation <- function(method = c("ssGSEA", "logmedian",
       labsize = labsize,
       titlesize = titlesize,
       pointSize = pointSize,
-      ignore_NAs = ignore_NAs
+      ignore_NAs = ignore_NAs,
+      p.adjust.method = p.adjust.method
     )
 
   } else if (method %in% c("ssGSEA", "logmedian", "ranking")) {
@@ -397,7 +409,8 @@ VariableAssociation <- function(method = c("ssGSEA", "logmedian",
       discrete_colors = discrete_colors,
       continuous_color = continuous_color,
       color_palette = color_palette,
-      printplt = printplt
+      printplt = printplt,
+      p.adjust.method = p.adjust.method
     )
   }
 
diff --git a/R/runGSEA.R b/R/runGSEA.R
index bed6e869..d80279c1 100644
--- a/R/runGSEA.R
+++ b/R/runGSEA.R
@@ -28,7 +28,13 @@
 #'   number of contrasts tested per signature and provides more stringent
 #'   control of false discovery rate across multiple comparisons. If `FALSE`,
 #'   the function only corrects for the number of gene sets.
-#'
+#'   
+#' @param p.adjust.method Character string specifying the method to use for
+#'   multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+#'   default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+#'   \code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+#'   Passed to \code{\link[stats]{p.adjust}}. 
+#'   
 #' @param nPermSimple Number of permutations in the simple fgsea implementation
 #'   for preliminary estimation of P-values. Parameter from fgsea.
 #'
@@ -57,7 +63,7 @@
 #'
 #' @importFrom fgsea fgsea
 #' @export
-runGSEA <- function(DEGList, gene_sets, stat = NULL, ContrastCorrection=FALSE, nPermSimple=10000) {
+runGSEA <- function(DEGList, gene_sets, stat = NULL, ContrastCorrection=FALSE, nPermSimple=10000, p.adjust.method="BH") {
 
   # Initialize storage for results across contrasts
   results_by_contrast <- list()
@@ -150,7 +156,7 @@ runGSEA <- function(DEGList, gene_sets, stat = NULL, ContrastCorrection=FALSE, n
     combined_df <- do.call(rbind, Map(cbind, results_by_contrast, df_name = names(results_by_contrast)))
 
     # Step 2: Adjust p-values across all data
-    combined_df$padj <- stats::p.adjust(combined_df$pval, method = "BH")
+    combined_df$padj <- stats::p.adjust(combined_df$pval,  method = p.adjust.method)
 
     # Step 3: Split back into the original list structure
     list_of_dfs <- split(combined_df, combined_df$df_name)
@@ -163,7 +169,7 @@ runGSEA <- function(DEGList, gene_sets, stat = NULL, ContrastCorrection=FALSE, n
 
     # Step 1: Adjust p-values for each data frame individually
     results_by_contrast <- lapply(results_by_contrast, function(df) {
-      df$padj <- stats::p.adjust(df$pval, method = "BH")  # Adjust p-values per data frame
+      df$padj <- stats::p.adjust(df$pval,  method = p.adjust.method)  # Adjust p-values per data frame
       return(df)
     })
 
diff --git a/README.Rmd b/README.Rmd
index 1a883557..94f25d9c 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -30,7 +30,7 @@ knitr::opts_chunk$set(
 
 > **To cite `markeR` please use:** 
 >
->   Martins-Silva R, Kaizeler A, Barbosa-Morais N (2025). _markeR: an R Toolkit for Evaluating Gene Sets as Phenotypic Markers_. Gulbenkian Institute for Molecular Medicine, Faculdade de Medicina, Universidade de Lisboa, Lisbon, Portugal. R package version 1.0.0, https://github.com/DiseaseTranscriptomicsLab/markeR.
+> Martins-Silva R, Kaizeler A, Barbosa-Morais NL (2025). _markeR: An R Toolkit for Evaluating Gene Signatures as Phenotypic Markers_. doi:10.18129/B9.bioc.markeR, R package version 1.1.2, https://bioconductor.org/packages/markeR.
  
 The folder `inst/Paper/` is in the **paper** branch and contains all scripts and materials used in the original `markeR` paper to reproduce analyses and figures.  You can browse it [here](https://github.com/DiseaseTranscriptomicsLab/markeR/tree/paper/inst/Paper).
 
@@ -52,6 +52,7 @@ The folder `inst/Paper/` is in the **paper** branch and contains all scripts and
   - [4. Visualisation and Evaluation](#4-visualisation-and-evaluation)  
   - [5. Individual Gene Exploration (Optional)](#5-individual-gene-exploration-optional)  
   - [6. Compare with Reference Gene Sets (Optional)](#6-compare-with-reference-gene-sets-optional)  
+- [Python Bridge](#python-bridge)
 - [Contact](#contact)
 
 
@@ -237,6 +238,16 @@ the ratio of the number of genes in common over the total number of genes in the
 
 Filters can be applied based on similarity thresholds (e.g., minimum Jaccard, OR, or Fisher's test p-value).
 
+## Python Bridge
+
+For users who prefer Python, a lightweight bridge is available in
+`python/` that allows calling any `markeR` function from a Python
+environment via [`rpy2`](https://rpy2.github.io/). It includes a tutorial
+workflow script and a generic command-line wrapper. See
+[`python/README.md`](inst/python/README.md) for installation
+instructions and usage examples.
+
+
 ## Contact 
 
 📩 For any questions or concerns, feel free to reach out:
diff --git a/README.md b/README.md
index a72fc8b3..0daabf81 100644
--- a/README.md
+++ b/README.md
@@ -24,11 +24,10 @@ across experimental and clinical phenotypes.
 
 > **To cite `markeR` please use:**
 >
-> Martins-Silva R, Kaizeler A, Barbosa-Morais N (2025). *markeR: an R
-> Toolkit for Evaluating Gene Sets as Phenotypic Markers*. Gulbenkian
-> Institute for Molecular Medicine, Faculdade de Medicina, Universidade
-> de Lisboa, Lisbon, Portugal. R package version 1.0.0,
-> <https://github.com/DiseaseTranscriptomicsLab/markeR>.
+> Martins-Silva R, Kaizeler A, Barbosa-Morais NL (2025). *markeR: An R
+> Toolkit for Evaluating Gene Signatures as Phenotypic Markers*.
+> <doi:10.18129/B9.bioc.markeR>, R package version 1.1.2,
+> <https://bioconductor.org/packages/markeR>.
 
 The folder `inst/Paper/` is in the **paper** branch and contains all
 scripts and materials used in the original `markeR` paper to reproduce
@@ -54,6 +53,7 @@ analyses and figures. You can browse it
     (Optional)](#5-individual-gene-exploration-optional)  
   - [6. Compare with Reference Gene Sets
     (Optional)](#6-compare-with-reference-gene-sets-optional)  
+- [Python Bridge](#python-bridge)
 - [Contact](#contact)
 
 ## Installation
@@ -285,6 +285,15 @@ metrics are implemented:
 Filters can be applied based on similarity thresholds (e.g., minimum
 Jaccard, OR, or Fisher’s test p-value).
 
+## Python Bridge
+
+For users who prefer Python, a lightweight bridge is available in
+`python/` that allows calling any `markeR` function from a Python
+environment via [`rpy2`](https://rpy2.github.io/). It includes a
+tutorial workflow script and a generic command-line wrapper. See
+[`python/README.md`](inst/python/README.md) for installation
+instructions and usage examples.
+
 ## Contact
 
 📩 For any questions or concerns, feel free to reach out:
diff --git a/man/CohenD_allConditions.Rd b/man/CohenD_allConditions.Rd
index 5c0ea77c..55871294 100644
--- a/man/CohenD_allConditions.Rd
+++ b/man/CohenD_allConditions.Rd
@@ -9,7 +9,8 @@ CohenD_allConditions(
   metadata,
   gene_sets,
   variable,
-  mode = c("simple", "medium", "extensive")
+  mode = c("simple", "medium", "extensive"),
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -37,6 +38,12 @@ groups.
 \item \code{"extensive"}: All possible groupwise contrasts, ensuring balance in the
 number of terms on each side.
 }}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A named list where each element corresponds to a gene signature. Each
diff --git a/man/CohenF_allConditions.Rd b/man/CohenF_allConditions.Rd
index 8a5f4b49..df25f8c1 100644
--- a/man/CohenF_allConditions.Rd
+++ b/man/CohenF_allConditions.Rd
@@ -4,7 +4,13 @@
 \alias{CohenF_allConditions}
 \title{Compute Cohen's f for All Gene Signatures Across a Categorical Variable}
 \usage{
-CohenF_allConditions(data, metadata, gene_sets, variable)
+CohenF_allConditions(
+  data,
+  metadata,
+  gene_sets,
+  variable,
+  p.adjust.method = "BH"
+)
 }
 \arguments{
 \item{data}{A data frame of gene expression data, with genes as rows and
@@ -21,6 +27,12 @@ downregulated).}
 
 \item{variable}{A string specifying the categorical variable in
 \code{metadata} used to model the gene signature scores.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A named list where each element corresponds to a gene signature. Each
diff --git a/man/FPR_Simulation.Rd b/man/FPR_Simulation.Rd
index b9fd5bf2..fb94bf7c 100644
--- a/man/FPR_Simulation.Rd
+++ b/man/FPR_Simulation.Rd
@@ -19,7 +19,8 @@ FPR_Simulation(
   mode = c("none", "simple", "medium", "extensive"),
   ColorValues = NULL,
   ncol = NULL,
-  nrow = NULL
+  nrow = NULL,
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -77,6 +78,12 @@ grid layout. If \code{NULL}, layout is auto-calculated.}
 
 \item{nrow}{Integer. Number of rows for arranging signature plots in a grid
 layout. If \code{NULL}, layout is auto-calculated.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 Invisibly returns a list containing:
diff --git a/man/GSEA_VariableAssociation.Rd b/man/GSEA_VariableAssociation.Rd
index f3015a47..915560d3 100644
--- a/man/GSEA_VariableAssociation.Rd
+++ b/man/GSEA_VariableAssociation.Rd
@@ -20,7 +20,8 @@ GSEA_VariableAssociation(
   titlesize = 14,
   pointSize = 5,
   ignore_NAs = FALSE,
-  printplt = TRUE
+  printplt = TRUE,
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -96,6 +97,12 @@ removed before analysis, leading to a loss of data to be fitted in the
 model.}
 
 \item{printplt}{Boolean specifying if plot is to be printed. Default: \code{TRUE}.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A list with two elements:
diff --git a/man/PlotScores.Rd b/man/PlotScores.Rd
index 8d17eb62..a685ddf2 100644
--- a/man/PlotScores.Rd
+++ b/man/PlotScores.Rd
@@ -31,7 +31,8 @@ PlotScores(
   sig_threshold = 0.05,
   cohen_threshold = 0.5,
   colorPalette = "Set3",
-  cor = c("pearson", "spearman", "kendall")
+  cor = c("pearson", "spearman", "kendall"),
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -169,6 +170,12 @@ is manually specified.}
 \item{cor}{Correlation method for numeric variables. One of \code{"pearson"}
 (default), \code{"spearman"}, or \code{"kendall"}. Only applies when the variable is
 numeric and \code{method != "all"}.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}. Only if \code{method == "all"}.}
 }
 \value{
 Depending on \code{method}:
diff --git a/man/Score_VariableAssociation.Rd b/man/Score_VariableAssociation.Rd
index 6c7c43ba..4e6a50df 100644
--- a/man/Score_VariableAssociation.Rd
+++ b/man/Score_VariableAssociation.Rd
@@ -23,7 +23,8 @@ Score_VariableAssociation(
   discrete_colors = NULL,
   continuous_color = "#8C6D03",
   color_palette = "Set2",
-  printplt = TRUE
+  printplt = TRUE,
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -85,6 +86,12 @@ variables. Default: \code{"#8C6D03"}.}
 variables. Default: \code{"Set2"}.}
 
 \item{printplt}{Boolean specifying if plot is to be printed. Default: \code{TRUE}.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A list with:
diff --git a/man/VariableAssociation.Rd b/man/VariableAssociation.Rd
index 6b22c5e6..8e055d58 100644
--- a/man/VariableAssociation.Rd
+++ b/man/VariableAssociation.Rd
@@ -24,7 +24,8 @@ VariableAssociation(
   discrete_colors = NULL,
   continuous_color = "#8C6D03",
   color_palette = "Set2",
-  printplt = TRUE
+  printplt = TRUE,
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -88,6 +89,12 @@ points (default: \code{"#8C6D03"}).}
 categorical variables (default: \code{"Set2"}).}
 
 \item{printplt}{Logical. If \code{TRUE}, plots are printed. Default: \code{TRUE}.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A list with method-specific results and ggplot2-based visualizations:
diff --git a/man/compute_stat_tests.Rd b/man/compute_stat_tests.Rd
index 9e2df313..8b290b58 100644
--- a/man/compute_stat_tests.Rd
+++ b/man/compute_stat_tests.Rd
@@ -10,7 +10,8 @@ compute_stat_tests(
   cols = NULL,
   numeric = "pearson",
   categorical_bin = "t.test",
-  categorical_multi = "anova"
+  categorical_multi = "anova",
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -30,6 +31,12 @@ Options: \code{"t.test"} (default) or \code{"wilcoxon"}.}
 \item{categorical_multi}{The statistical test for multi-level categorical
 variables.
 Options: \code{"anova"} (default) or \code{"kruskal-wallis"}.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A named list (one entry per variable being analysed) where each
diff --git a/man/runGSEA.Rd b/man/runGSEA.Rd
index c1df8ef7..a73336a4 100644
--- a/man/runGSEA.Rd
+++ b/man/runGSEA.Rd
@@ -9,7 +9,8 @@ runGSEA(
   gene_sets,
   stat = NULL,
   ContrastCorrection = FALSE,
-  nPermSimple = 10000
+  nPermSimple = 10000,
+  p.adjust.method = "BH"
 )
 }
 \arguments{
@@ -48,6 +49,12 @@ the function only corrects for the number of gene sets.}
 
 \item{nPermSimple}{Number of permutations in the simple fgsea implementation
 for preliminary estimation of P-values. Parameter from fgsea.}
+
+\item{p.adjust.method}{Character string specifying the method to use for
+multiple testing correction. Must be one of \code{"BH"} (Benjamini-Hochberg,
+default), \code{"holm"}, \code{"hommel"}, \code{"bonferroni"},
+\code{"BY"} (Benjamini-Yekutieli), \code{"fdr"}, or \code{"none"}.
+Passed to \code{\link[stats]{p.adjust}}.}
 }
 \value{
 A named list where each element corresponds to a contrast. Each
diff --git a/python/README.md b/python/README.md
new file mode 100644
index 00000000..a2f7b7b6
--- /dev/null
+++ b/python/README.md
@@ -0,0 +1,151 @@
+# markeR Python Bridge
+
+This workspace provides two simple Python helpers for using the
+Bioconductor R package **markeR** via `rpy2`.
+
+- `markeR_to_python.py` provides a minimal example workflow 
+  demonstrating how to load the example data and call a specific 
+  function, serving as a template for users who wish to structure 
+  their own analysis script in Python.
+- `run_marker_function.py` is a generic command‑line wrapper that can
+  invoke *any* function exported by the markeR package.
+
+## Prerequisites
+
+* R (>=4.5) installed and on your `PATH`.
+* A Python virtual environment.  A `requirements.txt` file is provided in
+  this folder listing the needed packages (`rpy2`, `pandas`, `numpy` plus
+  optional `ipython`/`jupyter` for notebook usage).
+  To set up the environment:
+
+  ```bash
+  python -m venv .venv
+  source .venv/bin/activate
+  pip install -r requirements.txt
+  ```
+
+  After activation you can run the helper scripts using `python` from the
+  same environment.
+
+## Quick start
+
+### 1. Run the tutorial workflow
+
+```bash
+python markeR_to_python.py --tutorial --output tutorial.png
+open tutorial.png
+```
+
+This executes an example from the markeR vignette; consult the
+original R tutorial (link below) for a step‑by‑step description of each
+analysis step.
+
+### 2. Call any markeR function
+
+```bash
+python run_marker_function.py PlotScores \
+  --data counts_example \
+  --metadata metadata_example \
+  --gene_sets genesets_example \
+  --Variable "Condition" \
+  --method logmedian \
+  --nrow 1 \
+  --output my_plot.png
+```
+
+All parameters are passed as `--name value` pairs.  Use `--help-function` or
+refer to the online reference manual for argument names.  The two scripts
+are described briefly below.
+
+---
+
+## Tutorial script (`markeR_to_python.py`)
+
+A minimal Python example illustrating how to call markeR from within a Python 
+environment. It loads the example datasets, runs a selected function, and 
+generates a corresponding plot. The script is intended as a template for users 
+who wish to structure and extend their own analysis workflows in Python.
+
+```bash
+python markeR_to_python.py --tutorial [--output file.png]
+```
+
+---
+
+## Flexible CLI caller (`run_marker_function.py`)
+
+This wrapper constructs and executes R code on the fly, so you can run any
+markeR function without _writing_ R.  It handles type conversion, data
+loading and optional PNG output.
+
+Basic syntax:
+
+```bash
+python run_marker_function.py FUNCTION_NAME [OPTIONS]
+```
+
+Useful options:
+
+* `--param value` – named argument for the R function
+* `--output filename.png` – capture plot output
+* `--width` / `--height` – PNG dimensions (pixels; default 800×600)
+* `--help-function` – show documentation link/usage hint
+* `--verbose` – display generated R code prior to execution
+
+Example (scores only, no plot):
+
+```bash
+python run_marker_function.py CalculateScores \
+  --data counts_example \
+  --metadata metadata_example \
+  --gene_sets genesets_example \
+  --method logmedian
+```
+
+Example saving a plot to PNG (with custom dimensions):
+
+```bash
+python run_marker_function.py PlotScores \
+  --data counts_example \
+  --metadata metadata_example \
+  --gene_sets genesets_example \
+  --Variable "Condition" \
+  --method logmedian \
+  --nrow 1 \
+  --width 800 \
+  --height 400 \
+  --output my_plot.png
+```
+---
+
+## Built‑in example data
+
+`counts_example`, `metadata_example` and `genesets_example` are loaded
+automatically and mirror the objects used in the tutorial.  They let you
+try commands without supplying your own datasets.
+
+---
+
+## Tips & troubleshooting
+
+* Add `--verbose` to see the exact R code being run – handy when a
+  parameter doesn’t behave as expected.
+* If `--output` fails, check write permissions and ensure the directory
+  exists.
+* Installation problems usually indicate R is missing; make sure `R` is
+  on your path before installing Python dependencies.
+
+---
+
+## References
+
+* [markeR on Bioconductor](https://bioconductor.org/packages/markeR)
+* [Official R tutorial](https://diseasetranscriptomicslab.github.io/markeR/articles/Article_BenchmarkingMode.html)
+* [Reference manual](https://diseasetranscriptomicslab.github.io/markeR/reference/)
+* Paper: https://www.biorxiv.org/content/10.64898/2025.12.05.692517
+
+---
+
+*Notes:* scripts install markeR via BiocManager if it’s not already present.
+Warnings from `ggplot2` (e.g. about `aes_string()`) are harmless and come
+from the package itself.
diff --git a/python/markeR_to_python.py b/python/markeR_to_python.py
new file mode 100644
index 00000000..a2e3b09f
--- /dev/null
+++ b/python/markeR_to_python.py
@@ -0,0 +1,321 @@
+"""Simple workflow for using the Bioconductor package `markeR` from Python via rpy2.
+
+source .venv/bin/activate
+pip install -r requirements.txt
+python markeR_to_python.py --tutorial   # now should execute successfully
+
+This script shows how to:
+
+1. Configure rpy2 and ensure an R environment is available.
+2. Install Bioconductor and markeR if not already installed.
+3. Load R functions into Python.
+4. Demonstrate a markeR analysis using example data.
+
+Notes:
+- You must have R (>=4.5) installed on your system.
+- Install the Python package `rpy2` in the same environment where this script runs:
+    pip install rpy2
+
+For more details on markeR see https://bioconductor.org/packages/markeR
+
+This module also includes utilities (`plot_r_expression`, `plot_r_function`)
+that open an R graphics device and capture plots as PNG files.  When run inside a
+Jupyter notebook the plots are automatically displayed inline; otherwise the
+images are saved to a temporary file whose path is printed.
+"""
+
+from __future__ import annotations
+
+import sys
+
+# check that required Python libraries are installed before proceeding
+_missing = []
+try:
+    import numpy  # used by rpy2 and examples
+except ImportError:  # pragma: no cover - dependency check
+    _missing.append("numpy")
+try:
+    import pandas  # examples use it for conversions
+except ImportError:  # pragma: no cover
+    _missing.append("pandas")
+try:
+    import rpy2  # primary bridge to R
+except ImportError:  # pragma: no cover
+    _missing.append("rpy2")
+
+if _missing:
+    sys.exit(
+        "The following Python packages are required but not installed: %s.\n"
+        "Please install them (e.g. `pip install -r requirements.txt`)." %
+        ", ".join(_missing)
+    )
+
+# rpy2 imports
+import rpy2.robjects as ro
+from rpy2.robjects import pandas2ri
+from rpy2.robjects import conversion
+from rpy2.robjects.packages import importr, isinstalled
+
+# note: pandas2ri.activate() is deprecated; we use conversion contexts when
+# converting.  helpers below wrap the recommended API.
+
+
+def _to_py(obj):
+    """Convert an R object to a pandas/numpy equivalent."""
+    with conversion.localconverter(ro.default_converter + pandas2ri.converter):
+        return conversion.rpy2py(obj)
+
+# utilities for inline plotting (e.g. in Jupyter notebooks)
+import os
+import tempfile
+try:
+    from IPython.display import Image, display
+    _HAS_IPYTHON = True
+except ImportError:  # not running in notebook
+    _HAS_IPYTHON = False
+
+
+def _r_open_png(width=800, height=600, filename=None):
+    """Start an R PNG device, returning the filename used."""
+    if filename is None:
+        filename = tempfile.mktemp(suffix=".png")
+    else:
+        # Convert to absolute path to ensure R saves to the intended location
+        filename = os.path.abspath(filename)
+    # ensure the path exists
+    dirname = os.path.dirname(filename)
+    if dirname and not os.path.isdir(dirname):
+        os.makedirs(dirname, exist_ok=True)
+    ro.r(f'png("{filename}", width={width}, height={height})')
+    return filename
+
+
+def _r_close_device():
+    """Close the active R graphics device."""
+    ro.r('dev.off()')
+
+
+def plot_r_expression(expr: str, width=800, height=600, filename=None, display_plot=True):
+    """Evaluate an R expression that produces a plot and optionally display it.
+
+    Parameters
+    ----------
+    expr : str
+        R code string that generates a plot when evaluated.
+    width, height : int
+        Dimensions for the PNG device in pixels.
+    filename : str or None
+        Path to which the image should be saved. If None a temporary file
+        will be created.
+    display_plot : bool
+        If True and running under IPython, display the resulting PNG
+        inline. Otherwise the path is printed.
+
+    Returns
+    -------
+    str
+        The path to the saved PNG file.
+    """
+    fname = _r_open_png(width=width, height=height, filename=filename)
+    ro.r(expr)
+    _r_close_device()
+    if display_plot and _HAS_IPYTHON:
+        display(Image(filename=fname))
+    else:
+        print(f"plot written to {fname}")
+    return fname
+
+
+def plot_r_function(func_name: str, *args, width=800, height=600, filename=None,
+                    display_plot=True, **kwargs):
+    """Call an R plotting function by name and save/display result.
+
+    Any positional and keyword arguments are converted to their R
+    equivalents by rpy2.
+
+    Example::
+
+        plot_r_function('PlotScores', data=counts_example, metadata=metadata_example,
+                        gene_sets=genesets_example)
+    """
+    fname = _r_open_png(width=width, height=height, filename=filename)
+    rfunc = ro.r[func_name]
+    # convert kwargs into ro objects (rpy2 handles this automatically)
+    rfunc(*args, **kwargs)
+    _r_close_device()
+    
+    # Try to display if requested; if in actual IPython/Jupyter context, display inline;
+    # otherwise just print the path.
+    if display_plot and _HAS_IPYTHON:
+        try:
+            # Check if we're actually in an interactive IPython shell (not just that it's installed)
+            from IPython import get_ipython
+            ipython = get_ipython()
+            if ipython is not None:
+                display(Image(filename=fname))
+            else:
+                print(f"plot written to {fname}")
+        except Exception:
+            # If anything goes wrong (file not ready, not in IPython context, etc.)
+            print(f"plot written to {fname}")
+    else:
+        print(f"plot written to {fname}")
+    return fname
+
+
+# ---------------------------------------------------------------------------
+# Utility functions
+# ---------------------------------------------------------------------------
+
+def ensure_bioc_installed() -> None:
+    """Install Bioconductor's package manager if it is not already present."""
+    biocinstaller = "BiocManager"
+    if not isinstalled(biocinstaller):
+        ro.r('install.packages("{0}")'.format(biocinstaller))
+    ro.r('suppressMessages(require({0}))'.format(biocinstaller))
+
+
+def install_markeR() -> None:
+    """Install the markeR package from Bioconductor if not already installed.
+
+    Uses BiocManager to perform the installation.  After running this function
+    the package should be loadable via `importr("markeR")`.
+    """
+    ensure_bioc_installed()
+    if not isinstalled("markeR"):
+        ro.r('BiocManager::install("markeR", ask=FALSE, update=FALSE)')
+    ro.r('library(markeR)')
+
+
+def get_markeR_functions() -> ro.Environment:
+    """Return the markeR namespace so that functions can be accessed conveniently.
+
+    Example:
+        mark = get_markeR_functions()
+        scores = mark.CalculateScores(data=counts, metadata=metadata, gene_sets=genesets, method="logmedian")
+    """
+    install_markeR()
+    # importing via importr is more reliable than accessing `ro.r['markeR']`.
+    try:
+        return importr("markeR")
+    except Exception as e:
+        raise RuntimeError("Unable to load markeR package: %s" % e)
+
+
+
+
+
+# ---------------------------------------------------------------------------
+# Tutorial helpers using markeR example data
+# ---------------------------------------------------------------------------
+
+def load_benchmark_examples():
+    """Load the built-in example data and gene sets from the markeR package.
+
+    Returns a tuple `(counts, metadata, genesets)` where each element is an R
+    object.  You can convert them to pandas objects if desired.
+    """
+    # ensure package is installed and loaded
+    install_markeR()
+    # load the three example datasets provided by the vignette
+    ro.r('data("genesets_example", package="markeR")')
+    ro.r('data("counts_example", package="markeR")')
+    ro.r('data("metadata_example", package="markeR")')
+    genesets = ro.r('genesets_example')
+    counts = ro.r('counts_example')
+    metadata = ro.r('metadata_example')
+    return counts, metadata, genesets
+
+
+def tutorial_benchmark(output_file=None):
+    """Demonstrate a small benchmarking mode example from the markeR vignette.
+    
+    Parameters
+    ----------
+    output_file : str or None
+        If provided, saves the display output to a file using the R graphics device.
+        Note: The markeR::PlotScores function outputs to the active graphics device.
+    """
+    print("-- loading example data from markeR")
+    counts, metadata, genesets = load_benchmark_examples()
+
+    # show dimensions of the data
+    print("counts matrix dimensions:", ro.r('dim')(counts))
+    print("metadata dimensions:", ro.r('dim')(metadata))
+    print("available gene sets:", list(genesets.names))
+
+    # run CalculateScores (logmedian method) as in the tutorial
+    calculate = ro.r['CalculateScores']
+    print("-- calculating scores using logmedian method")
+    df_scores = calculate(data=counts,
+                          metadata=metadata,
+                          method="logmedian",
+                          gene_sets=genesets)
+
+    # df_scores is an R list with one element per gene set; convert first one
+    # to pandas for display
+    first_name = list(df_scores.names)[0]
+    r_first = df_scores.rx2(first_name)
+    try:
+        import pandas as pd
+        pd_first = _to_py(r_first)
+        print(f"first gene set ({first_name}) scores (head):\n", pd_first.head())
+    except ImportError:
+        print("pandas not available; skipping conversion of results to DataFrame")
+
+    # Generate plot using PlotScores
+    # Note: markeR's PlotScores function creates an interactive plot or writes to the current device
+    print("-- generating a simple score plot")
+    
+    if output_file:
+        # Set up PNG device
+        output_file = os.path.abspath(output_file)
+        dirname = os.path.dirname(output_file)
+        if dirname and not os.path.isdir(dirname):
+            os.makedirs(dirname, exist_ok=True)
+        ro.r(f'png("{output_file}", width=800, height=400)')
+        print(f"   saving to: {output_file}")
+    
+    # Call the plotting function and force evaluation of the returned plot
+    ro.r('''
+    p <- PlotScores(
+        data = counts_example,
+        metadata = metadata_example,
+        gene_sets = genesets_example,
+        Variable = "Condition",
+        method = "logmedian",
+        nrow=1
+    )
+    print(p)
+    ''')
+    
+    if output_file:
+        # Close device
+        ro.r('dev.off()')
+        print(f"plot saved to {output_file}")
+
+# ---------------------------------------------------------------------------
+# Command-line interface
+# ---------------------------------------------------------------------------
+
+if __name__ == "__main__":
+    output_file = None
+    
+    # Parse command line arguments
+    if len(sys.argv) > 1:
+        if sys.argv[1] == "--tutorial":
+            # Check if --output flag is present
+            if len(sys.argv) > 2 and sys.argv[2].startswith("--output"):
+                if sys.argv[2] == "--output" and len(sys.argv) > 3:
+                    output_file = sys.argv[3]
+                elif "=" in sys.argv[2]:
+                    output_file = sys.argv[2].split("=", 1)[1]
+            tutorial_benchmark(output_file=output_file)
+        else:
+            print("usage: python markeR_to_python.py --tutorial [--output FILENAME]")
+            print("  --tutorial              : load markeR example data and compute logmedian scores")
+            print("  --output FILENAME       : save plot to specified PNG file (optional)")
+            print("                            example: python markeR_to_python.py --tutorial --output my_plot.png")
+    else:
+        print("usage: python markeR_to_python.py --tutorial [--output FILENAME]")
+        print("See the module docstring for more details.")
\ No newline at end of file
diff --git a/python/requirements.txt b/python/requirements.txt
new file mode 100644
index 00000000..fdb0e44a
--- /dev/null
+++ b/python/requirements.txt
@@ -0,0 +1,9 @@
+# install packages for markeR_to_python script:
+#   pip install -r requirements.txt
+
+rpy2>=3.6
+pandas
+numpy
+# optional (for notebook inline display)
+ipython
+jupyter
diff --git a/python/run_marker_function.py b/python/run_marker_function.py
new file mode 100644
index 00000000..5ca33375
--- /dev/null
+++ b/python/run_marker_function.py
@@ -0,0 +1,366 @@
+"""Flexible wrapper to call any markeR R function from Python.
+
+This script allows you to call markeR functions directly without writing Python code.
+
+Usage examples:
+    python run_marker_function.py PlotScores --help-function
+    python run_marker_function.py CalculateScores \\
+        --data counts_example --metadata metadata_example \\
+        --gene_sets genesets_example --method logmedian --verbose
+    
+    python run_marker_function.py PlotScores \\
+        --data counts_example --metadata metadata_example \\
+        --gene_sets genesets_example --Variable "Condition" \\
+        --method logmedian --nrow 1 --output plot.png
+
+Options:
+    --help-function      Show R documentation for the function
+    --verbose            Print the generated R code before executing
+    --output FILE        Save plot to PNG file
+
+For built-in example data, use the names: counts_example, metadata_example, genesets_example
+"""
+
+import sys
+import argparse
+import json
+import os
+import re
+
+# Check dependencies
+_missing = []
+try:
+    import rpy2.robjects as ro
+    from rpy2.robjects.packages import importr, isinstalled
+except ImportError:
+    _missing.append("rpy2")
+
+if _missing:
+    sys.exit(
+        "The following Python packages are required but not installed: %s.\n"
+        "Please install them (e.g. `pip install rpy2`)." % ", ".join(_missing)
+    )
+
+def ensure_bioc_installed() -> None:
+    """Install Bioconductor's package manager if it is not already present."""
+    biocinstaller = "BiocManager"
+    if not isinstalled(biocinstaller):
+        ro.r('install.packages("{0}")'.format(biocinstaller))
+    ro.r('suppressMessages(require({0}))'.format(biocinstaller))
+
+
+def install_markeR() -> None:
+    """Install the markeR package from Bioconductor if not already installed."""
+    ensure_bioc_installed()
+    if not isinstalled("markeR"):
+        ro.r('BiocManager::install("markeR", ask=FALSE, update=FALSE)')
+    ro.r('library(markeR)')
+
+
+def load_example_data():
+    """Load built-in markeR example datasets into R namespace."""
+    install_markeR()
+    # Load the example datasets
+    ro.r('data("genesets_example", package="markeR")')
+    ro.r('data("counts_example", package="markeR")')
+    ro.r('data("metadata_example", package="markeR")')
+    print("Loaded markeR example datasets: counts_example, metadata_example, genesets_example")
+
+
+def parse_parameter(value: str):
+    """
+    Parse a parameter value intelligently.
+    - Numbers become numeric
+    - "true"/"false" become logical
+    - "null" becomes NULL
+    - R object names (e.g., counts_example) are kept as-is
+    - JSON objects/arrays become R equivalents
+    - Strings are kept as strings
+    """
+    value_lower = value.lower()
+    
+    # Handle boolean
+    if value_lower == "true":
+        return "TRUE"
+    if value_lower == "false":
+        return "FALSE"
+    if value_lower == "null":
+        return "NULL"
+    
+    # Handle numbers
+    try:
+        if "." in value:
+            float(value)
+            return value
+        else:
+            int(value)
+            return value
+    except ValueError:
+        pass
+    
+    # Check if it's a known R object name (example data)
+    known_objects = ["counts_example", "metadata_example", "genesets_example"]
+    if value in known_objects:
+        return value
+    
+    # Try JSON parsing for objects/arrays
+    try:
+        json.loads(value)
+        # If it parses as JSON, return as-is (user can provide lists as JSON)
+        return value
+    except (json.JSONDecodeError, ValueError):
+        pass
+    
+    # Default: treat as string, escaping any internal quotes
+    value = value.replace('"', '\\"')
+    return f'"{value}"'
+
+
+def build_r_call(function_name: str, params: dict, output_file: str = None, width: int = 800, height: int = 600) -> str:
+    """
+    Build an R function call string from parameters.
+    
+    Parameters
+    ----------
+    function_name : str
+        Name of the R function to call
+    params : dict
+        Dictionary of parameter names and values
+    output_file : str
+        If provided, set up PNG device before the call and close after
+    width : int
+        PNG width in pixels (default: 800)
+    height : int
+        PNG height in pixels (default: 600)
+        
+    Returns
+    -------
+    str
+        Complete R code to execute
+    """
+    # Remove output_file from params if present
+    params = {k: v for k, v in params.items() if k != "output_file"}
+    
+    # Build parameter list
+    param_strings = []
+    for key, value in params.items():
+        parsed_value = parse_parameter(value)
+        param_strings.append(f"{key} = {parsed_value}")
+    
+    param_str = ", ".join(param_strings)
+    
+    # Build R code
+    r_code = ""
+    
+    if output_file:
+        output_file = os.path.abspath(output_file)
+        dirname = os.path.dirname(output_file)
+        if dirname and not os.path.isdir(dirname):
+            os.makedirs(dirname, exist_ok=True)
+        r_code += f'png("{output_file}", width={width}, height={height})\n'
+        r_code += f"result <- {function_name}({param_str})\n"
+        r_code += "tryCatch(print(result), error=function(e) { invisible(NULL) })\n"
+        r_code += "dev.off()\n"
+    else:
+        # For screen output, try to print the result
+        r_code += f"result <- {function_name}({param_str})\n"
+        r_code += "tryCatch(print(result), error=function(e) { cat('Function executed.\\n') })\n"
+    
+    return r_code
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Call any markeR R function from Python",
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog="""
+Examples:
+  python run_marker_function.py CalculateScores \\
+    --data counts_example --metadata metadata_example \\
+    --gene_sets genesets_example --method logmedian
+    
+  python run_marker_function.py PlotScores \\
+    --data counts_example --metadata metadata_example \\
+    --gene_sets genesets_example --Variable "Condition" \\
+    --method logmedian --output my_plot.png
+
+Built-in example data: counts_example, metadata_example, genesets_example
+        """
+    )
+    
+    parser.add_argument(
+        "function_name",
+        help="Name of the markeR function to call (e.g., CalculateScores, PlotScores)"
+    )
+    
+    parser.add_argument(
+        "--help-function",
+        action="store_true",
+        help="Show help for the R function (instead of calling it)"
+    )
+    
+    parser.add_argument(
+        "--verbose",
+        action="store_true",
+        help="Print the generated R code before executing"
+    )
+    
+    parser.add_argument(
+        "--output",
+        help="Save plot output to a PNG file"
+    )
+    
+    parser.add_argument(
+        "--width",
+        type=int,
+        default=800,
+        help="PNG width in pixels (default: 800)"
+    )
+    
+    parser.add_argument(
+        "--height",
+        type=int,
+        default=600,
+        help="PNG height in pixels (default: 600)"
+    )
+    
+    # Allow arbitrary parameters
+    parser.add_argument(
+        "params",
+        nargs="*",
+        help="Parameters as --name value pairs (e.g., --data counts_example --method logmedian)"
+    )
+    
+    # Handle --help for specific functions
+    if len(sys.argv) > 1 and sys.argv[1] not in ["--help", "-h"]:
+        if "--help-function" in sys.argv:
+            func_name = sys.argv[1]
+            print(f"\n{'='*70}")
+            print(f"Help for markeR::{func_name}")
+            print(f"{'='*70}\n")
+            install_markeR()
+            # Try to display help
+            try:
+                # Get function signature and description
+                ro.r(f'''
+library(markeR)
+cat("Function: {func_name}\\n\\n")
+# Try to get help
+tryCatch({{
+  help_file <- help("{func_name}", package="markeR")
+  # Get description from help
+}}, error = function(e) {{
+  cat("Help available at: https://diseasetranscriptomicslab.github.io/markeR/reference/{func_name}.html\\n")
+}})
+''')
+            except Exception as e:
+                pass
+            
+            print(f"\nDocumentation:")
+            print(f"  https://diseasetranscriptomicslab.github.io/markeR/reference/{func_name}.html")
+            print(f"\nTo use this function:")
+            print(f"  python run_marker_function.py {func_name} --param1 value1 --param2 value2 [--output output.png]")
+            print(f"\nTip: Use --verbose flag to see generated R code")
+            print(f"  python run_marker_function.py {func_name} --verbose --param1 value1 ...\n")
+            print(f"{'='*70}\n")
+            return
+    
+    # Parse args
+    if len(sys.argv) < 2:
+        parser.print_help()
+        return
+    
+    func_name = sys.argv[1]
+    
+    # Validate function name to prevent code injection
+    if not re.match(r'^[A-Za-z][A-Za-z0-9_.]*$', func_name):
+        sys.exit(f"Error: invalid function name '{func_name}'. "
+                 "Function names must start with a letter and contain only letters, digits, dots or underscores.")
+    
+    output_file = None
+    width = 800
+    height = 600
+    
+    # Parse remaining arguments as key-value pairs
+    params = {}
+    i = 2
+    while i < len(sys.argv):
+        arg = sys.argv[i]
+        if arg == "--output" and i + 1 < len(sys.argv):
+            output_file = sys.argv[i + 1]
+            i += 2
+        elif arg == "--width" and i + 1 < len(sys.argv):
+            try:
+                width = int(sys.argv[i + 1])
+            except ValueError:
+                print(f"Error: --width must be a number, got '{sys.argv[i + 1]}'")
+                sys.exit(1)
+            i += 2
+        elif arg == "--height" and i + 1 < len(sys.argv):
+            try:
+                height = int(sys.argv[i + 1])
+            except ValueError:
+                print(f"Error: --height must be a number, got '{sys.argv[i + 1]}'")
+                sys.exit(1)
+            i += 2
+        elif arg in ["--verbose", "--help-function"]:
+            # Skip flags that are not parameters
+            i += 1
+        elif arg.startswith("--"):
+            key = arg[2:]  # Remove --
+            if i + 1 < len(sys.argv) and not sys.argv[i + 1].startswith("--"):
+                value = sys.argv[i + 1]
+                params[key] = value
+                i += 2
+            else:
+                # Boolean flag
+                params[key] = "TRUE"
+                i += 1
+        else:
+            i += 1
+    
+    # Ensure markeR is installed
+    print("Installing markeR if needed...")
+    install_markeR()
+    load_example_data()
+    
+    # Build and execute the R call
+    print(f"\nCalling {func_name} with parameters:")
+    for key, value in params.items():
+        print(f"  {key} = {value}")
+    
+    if output_file:
+        print(f"  Saving plot to: {output_file}")
+        print(f"  PNG dimensions: {width}x{height} pixels")
+    
+    r_code = build_r_call(func_name, params, output_file, width, height)
+    
+    # Show R code if verbose mode
+    if "--verbose" in sys.argv:
+        print(f"\n{'='*70}")
+        print("Generated R code:")
+        print(f"{'='*70}")
+        print(r_code)
+        print(f"{'='*70}\n")
+    
+    print(f"Executing R code...\n")
+    print("=" * 60)
+    
+    try:
+        ro.r(r_code)
+        print("=" * 60)
+        if output_file:
+            print(f"\n✓ Plot saved to: {output_file}")
+        else:
+            print(f"\n✓ Function executed successfully")
+    except Exception as e:
+        print("=" * 60)
+        print(f"\n✗ Error executing function: {e}")
+        if "--verbose" not in sys.argv:
+            print("\nTip: Use --verbose flag to see the generated R code")
+            print(f"  python run_marker_function.py {func_name} --verbose [other options]")
+        sys.exit(1)
+
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
new file mode 100644
index 00000000..49c490f4
Binary files /dev/null and b/tests/testthat/Rplots.pdf differ