Tighten contrast facades and add lmer support.

Require the aggregate-only facades to run on protein-level LFQData, add a peptide-level lmer facade, and align the vignettes and pkgdown metadata with the stricter facade contracts.

Made-with: Cursor

Author: wolski

DESCRIPTION

@@ -23,7 +23,8 @@ Encoding: UTF-8
 LazyData: true
 LazyDataCompression: xz
 VignetteBuilder: knitr
-URL: https://github.com/wolski/prolfqua
+URL: https://github.com/wolski/prolfqua,
+    https://wolski.github.io/prolfqua/
 BugReports: https://github.com/wolski/prolfqua/issues
 Imports:
     bookdown,

NAMESPACE

@@ -12,6 +12,7 @@ export(ContrastsLMFacade)
 export(ContrastsLMMissingFacade)
 export(ContrastsLimma)
 export(ContrastsLimmaFacade)
+export(ContrastsLmerFacade)
 export(ContrastsMissing)
 export(ContrastsModerated)
 export(ContrastsModeratedDEqMS)

R/ContrastsFacades.R

@@ -1,5 +1,40 @@
 # ContrastsFacades -----
 
+.assert_aggregated_facade_input <- function(lfqdata, facade_name) {
+  subject_id <- lfqdata$subject_Id()
+  hierarchy_keys <- lfqdata$config$hierarchy_keys()
+  if (!identical(subject_id, hierarchy_keys)) {
+    stop(
+      facade_name,
+      " requires aggregated LFQData. ",
+      "`lfqdata$subject_Id()` must equal `lfqdata$config$hierarchy_keys()`. ",
+      "Aggregate first.",
+      call. = FALSE
+    )
+  }
+}
+
+.assert_nested_facade_input <- function(lfqdata, facade_name) {
+  subject_id <- lfqdata$subject_Id()
+  hierarchy_keys <- lfqdata$config$hierarchy_keys()
+  if (!(all(subject_id %in% hierarchy_keys) && length(subject_id) < length(hierarchy_keys))) {
+    stop(
+      facade_name,
+      " requires LFQData with additional hierarchy below `subject_Id()`. ",
+      "`lfqdata$subject_Id()` must be a strict subset of ",
+      "`lfqdata$config$hierarchy_keys()`. Do not aggregate first.",
+      call. = FALSE
+    )
+  }
+}
+
+.add_facade_column <- function(res, facade_name) {
+  if (!("facade" %in% colnames(res))) {
+    res <- dplyr::mutate(res, facade = facade_name, .before = 1)
+  }
+  res
+}
+
 #' Limma contrast analysis facade
 #'
 #' Encapsulates the pipeline: \code{\link{strategy_limma}} ->
@@ -8,9 +43,9 @@
 #' @export
 #' @family modelling
 #' @examples
-#' istar <- sim_lfq_data_peptide_config()
+#' istar <- sim_lfq_data_protein_config()
 #' lfqdata <- LFQData$new(istar$data, istar$config)
-#' lfqdata <- lfqdata$get_Transformer()$log2()$lfq
+#' lfqdata$rename_response("transformedIntensity")
 #' contrasts <- c("A_vs_Ctrl" = "group_A - group_Ctrl")
 #' fa <- ContrastsLimmaFacade$new(lfqdata, "~ group_", contrasts)
 #' head(fa$get_contrasts())
@@ -29,6 +64,7 @@ ContrastsLimmaFacade <- R6::R6Class(
     #' @param contrasts named character vector of contrasts
     #' @param ... passed to \code{\link{strategy_limma}} (e.g. trend, robust)
     initialize = function(lfqdata, modelstr, contrasts, ...) {
+      .assert_aggregated_facade_input(lfqdata, "ContrastsLimmaFacade")
       response <- lfqdata$config$get_response()
       full_formula <- paste(response, modelstr)
       strat          <- strategy_limma(full_formula, ...)
@@ -37,7 +73,9 @@ ContrastsLimmaFacade <- R6::R6Class(
     },
     #' @description get contrast results
     #' @param ... passed to ContrastsLimma$get_contrasts
-    get_contrasts = function(...) self$contrast$get_contrasts(...),
+    get_contrasts = function(...) {
+      .add_facade_column(self$contrast$get_contrasts(...), "limma")
+    },
     #' @description get ContrastsPlotter
     #' @param ... passed to ContrastsLimma$get_Plotter
     get_Plotter   = function(...) self$contrast$get_Plotter(...),
@@ -57,9 +95,9 @@ ContrastsLimmaFacade <- R6::R6Class(
 #' @export
 #' @family modelling
 #' @examples
-#' istar <- sim_lfq_data_peptide_config()
+#' istar <- sim_lfq_data_protein_config()
 #' lfqdata <- LFQData$new(istar$data, istar$config)
-#' lfqdata <- lfqdata$get_Transformer()$log2()$lfq
+#' lfqdata$rename_response("transformedIntensity")
 #' contrasts <- c("A_vs_Ctrl" = "group_A - group_Ctrl")
 #' fa <- ContrastsLMFacade$new(lfqdata, "~ group_", contrasts)
 #' head(fa$get_contrasts())
@@ -78,6 +116,7 @@ ContrastsLMFacade <- R6::R6Class(
     #' @param contrasts named character vector of contrasts
     #' @param ... passed to \code{\link{strategy_lm}}
     initialize = function(lfqdata, modelstr, contrasts, ...) {
+      .assert_aggregated_facade_input(lfqdata, "ContrastsLMFacade")
       response <- lfqdata$config$get_response()
       full_formula <- paste(response, modelstr)
       strat          <- strategy_lm(full_formula, ...)
@@ -86,7 +125,9 @@ ContrastsLMFacade <- R6::R6Class(
     },
     #' @description get contrast results
     #' @param ... passed to ContrastsModerated$get_contrasts
-    get_contrasts = function(...) self$contrast$get_contrasts(...),
+    get_contrasts = function(...) {
+      .add_facade_column(self$contrast$get_contrasts(...), "lm")
+    },
     #' @description get ContrastsPlotter
     #' @param ... passed to ContrastsModerated$get_Plotter
     get_Plotter   = function(...) self$contrast$get_Plotter(...),
@@ -97,6 +138,66 @@ ContrastsLMFacade <- R6::R6Class(
 )
 
 
+#' Lmer contrast analysis facade
+#'
+#' Encapsulates the pipeline: \code{\link{strategy_lmer}} ->
+#' \code{\link{build_model}} -> \code{\link{Contrasts}} ->
+#' \code{\link{ContrastsModerated}}.
+#'
+#' This facade requires data with hierarchy below the analysis subject, for
+#' example peptide-level measurements nested within proteins.
+#'
+#' @export
+#' @family modelling
+#' @examples
+#' istar <- sim_lfq_data_peptide_config()
+#' istar$config <- old2new(istar$config)
+#' lfqdata <- LFQData$new(istar$data, istar$config)
+#' lfqdata <- lfqdata$get_Transformer()$log2()$lfq
+#' contrasts <- c("A_vs_Ctrl" = "group_A - group_Ctrl")
+#' fa <- ContrastsLmerFacade$new(
+#'   lfqdata,
+#'   "~ group_ + (1 | peptide_Id) + (1 | sampleName)",
+#'   contrasts
+#' )
+#' head(fa$get_contrasts())
+#' fa$to_wide()
+ContrastsLmerFacade <- R6::R6Class(
+  "ContrastsLmerFacade",
+  public = list(
+    #' @field model Model object
+    model = NULL,
+    #' @field contrast ContrastsModerated object
+    contrast = NULL,
+    #' @description
+    #' initialize
+    #' @param lfqdata LFQData object
+    #' @param modelstr model formula string (e.g. "~ group_ + (1 | peptide_Id)")
+    #' @param contrasts named character vector of contrasts
+    #' @param ... passed to \code{\link{strategy_lmer}}
+    initialize = function(lfqdata, modelstr, contrasts, ...) {
+      .assert_nested_facade_input(lfqdata, "ContrastsLmerFacade")
+      response <- lfqdata$config$get_response()
+      full_formula <- paste(response, modelstr)
+      strat          <- strategy_lmer(full_formula, ...)
+      self$model     <- build_model(lfqdata, strat)
+      self$contrast  <- ContrastsModerated$new(Contrasts$new(self$model, contrasts))
+    },
+    #' @description get contrast results
+    #' @param ... passed to ContrastsModerated$get_contrasts
+    get_contrasts = function(...) {
+      .add_facade_column(self$contrast$get_contrasts(...), "lmer")
+    },
+    #' @description get ContrastsPlotter
+    #' @param ... passed to ContrastsModerated$get_Plotter
+    get_Plotter   = function(...) self$contrast$get_Plotter(...),
+    #' @description convert results to wide format
+    #' @param ... passed to ContrastsModerated$to_wide
+    to_wide       = function(...) self$contrast$to_wide(...)
+  )
+)
+
+#'
 #' LM + missing-value imputation contrast analysis facade
 #'
 #' Encapsulates the pipeline: \code{\link{strategy_lm}} ->
@@ -136,6 +237,7 @@ ContrastsLMMissingFacade <- R6::R6Class(
     #' @param contrasts named character vector of contrasts
     #' @param ... passed to \code{\link{strategy_lm}}
     initialize = function(lfqdata, modelstr, contrasts, ...) {
+      .assert_aggregated_facade_input(lfqdata, "ContrastsLMMissingFacade")
       response <- lfqdata$config$get_response()
       full_formula <- paste(response, modelstr)
       strat                  <- strategy_lm(full_formula, ...)
@@ -147,7 +249,9 @@ ContrastsLMMissingFacade <- R6::R6Class(
     },
     #' @description get contrast results
     #' @param ... passed to ContrastsTable$get_contrasts
-    get_contrasts = function(...) self$contrast$get_contrasts(...),
+    get_contrasts = function(...) {
+      .add_facade_column(self$contrast$get_contrasts(...), "lm_missing")
+    },
     #' @description get ContrastsPlotter
     #' @param ... passed to ContrastsTable$get_Plotter
     get_Plotter   = function(...) self$contrast$get_Plotter(...),
@@ -189,6 +293,7 @@ ContrastsDEqMSFacade <- R6::R6Class(
     #' @param contrasts named character vector of contrasts
     #' @param ... passed to \code{\link{strategy_lm}}
     initialize = function(lfqdata, modelstr, contrasts, ...) {
+      .assert_aggregated_facade_input(lfqdata, "ContrastsDEqMSFacade")
       response <- lfqdata$config$get_response()
       full_formula <- paste(response, modelstr)
       strat         <- strategy_lm(full_formula, ...)
@@ -203,7 +308,9 @@ ContrastsDEqMSFacade <- R6::R6Class(
     },
     #' @description get contrast results
     #' @param ... passed to ContrastsModeratedDEqMS$get_contrasts
-    get_contrasts = function(...) self$contrast$get_contrasts(...),
+    get_contrasts = function(...) {
+      .add_facade_column(self$contrast$get_contrasts(...), "deqms")
+    },
     #' @description get ContrastsPlotter
     #' @param ... passed to ContrastsModeratedDEqMS$get_Plotter
     get_Plotter   = function(...) self$contrast$get_Plotter(...),
@@ -249,6 +356,7 @@ ContrastsROPECAFacade <- R6::R6Class(
     #' @param contrasts named character vector of contrasts
     #' @param ... passed to \code{\link{strategy_lm}}
     initialize = function(lfqdata, modelstr, contrasts, ...) {
+      .assert_nested_facade_input(lfqdata, "ContrastsROPECAFacade")
       response <- lfqdata$config$get_response()
       full_formula <- paste(response, modelstr)
       strat         <- strategy_lm(full_formula, ...)
@@ -293,7 +401,8 @@ ContrastsROPECAFacade <- R6::R6Class(
       standard_cols <- c(protein_Id, "modelName", "contrast", "avgAbd",
                          "diff", "FDR", "statistic", "std.error", "df",
                          "p.value", "conf.low", "conf.high", "sigma")
-      res[, standard_cols, drop = FALSE]
+      res <- res[, standard_cols, drop = FALSE]
+      .add_facade_column(res, "ropeca")
     },
     #' @description get ContrastsPlotter (uses standardized column names)
     #' @param FCthreshold fold change threshold

R/build_contrast_analysis.R

@@ -12,19 +12,20 @@
 #'   \code{lfqdata$config$get_response()}.
 #' @param contrasts named character vector of contrasts
 #'   (e.g. \code{c("A_vs_B" = "group_A - group_B")})
-#' @param method one of \code{"lm"}, \code{"lm_missing"}, \code{"limma"},
-#'   \code{"deqms"}, \code{"ropeca"}
+#' @param method one of \code{"lm"}, \code{"lmer"}, \code{"lm_missing"},
+#'   \code{"limma"}, \code{"deqms"}, \code{"ropeca"}
 #' @param ... additional arguments forwarded to the underlying strategy function
 #'   (e.g. \code{trend}, \code{robust} for \code{strategy_limma})
 #' @return one of \code{\link{ContrastsLimmaFacade}},
-#'   \code{\link{ContrastsLMFacade}}, \code{\link{ContrastsLMMissingFacade}},
-#'   \code{\link{ContrastsDEqMSFacade}}, or \code{\link{ContrastsROPECAFacade}}
+#'   \code{\link{ContrastsLMFacade}}, \code{\link{ContrastsLmerFacade}},
+#'   \code{\link{ContrastsLMMissingFacade}}, \code{\link{ContrastsDEqMSFacade}},
+#'   or \code{\link{ContrastsROPECAFacade}}
 #' @export
 #' @family modelling
 #' @examples
-#' istar <- sim_lfq_data_peptide_config()
+#' istar <- sim_lfq_data_protein_config(Nprot = 20)
 #' lfqdata <- LFQData$new(istar$data, istar$config)
-#' lfqdata <- lfqdata$get_Transformer()$log2()$lfq
+#' lfqdata$rename_response("transformedIntensity")
 #' contrasts <- c("A_vs_Ctrl" = "group_A - group_Ctrl")
 #'
 #' fa_lm    <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "lm")
@@ -33,26 +34,36 @@
 #' fa_limma <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "limma")
 #' head(fa_limma$get_contrasts())
 #'
-#' fa_ropeca <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "ropeca")
-#' head(fa_ropeca$get_contrasts())
-#'
-#' # lm_missing requires protein-level data
-#' istar_prot <- sim_lfq_data_protein_config(Nprot = 20)
-#' lfqdata_prot <- LFQData$new(istar_prot$data, istar_prot$config)
-#' lfqdata_prot$rename_response("transformedIntensity")
-#' fa_miss <- build_contrast_analysis(lfqdata_prot, "~ group_", contrasts, method = "lm_missing")
+#' fa_miss <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "lm_missing")
 #' head(fa_miss$get_contrasts())
 #'
-#' fa_deqms <- build_contrast_analysis(lfqdata_prot, "~ group_", contrasts, method = "deqms")
+#' fa_deqms <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "deqms")
 #' head(fa_deqms$get_contrasts())
+#'
+#' istar_pep <- sim_lfq_data_peptide_config()
+#' istar_pep$config <- old2new(istar_pep$config)
+#' lfqdata_pep <- LFQData$new(istar_pep$data, istar_pep$config)
+#' lfqdata_pep <- lfqdata_pep$get_Transformer()$log2()$lfq
+#'
+#' fa_lmer <- build_contrast_analysis(
+#'   lfqdata_pep,
+#'   "~ group_ + (1 | peptide_Id) + (1 | sampleName)",
+#'   contrasts,
+#'   method = "lmer"
+#' )
+#' head(fa_lmer$get_contrasts())
+#'
+#' fa_ropeca <- build_contrast_analysis(lfqdata_pep, "~ group_", contrasts, method = "ropeca")
+#' head(fa_ropeca$get_contrasts())
 build_contrast_analysis <- function(lfqdata,
                                     modelstr,
                                     contrasts,
-                                    method = c("lm", "lm_missing", "limma", "deqms", "ropeca"),
+                                    method = c("lm", "lmer", "lm_missing", "limma", "deqms", "ropeca"),
                                     ...) {
   method <- match.arg(method)
   switch(method,
     lm         = ContrastsLMFacade$new(lfqdata, modelstr, contrasts, ...),
+    lmer       = ContrastsLmerFacade$new(lfqdata, modelstr, contrasts, ...),
     lm_missing = ContrastsLMMissingFacade$new(lfqdata, modelstr, contrasts, ...),
     limma      = ContrastsLimmaFacade$new(lfqdata, modelstr, contrasts, ...),
     deqms      = ContrastsDEqMSFacade$new(lfqdata, modelstr, contrasts, ...),

README.md

@@ -3,7 +3,7 @@
 [![DOI](https://img.shields.io/badge/DOI-10.1021%2Facs.jproteome.2c00441-blue)](https://pubs.acs.org/doi/10.1021/acs.jproteome.2c00441)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.18775726.svg)](https://doi.org/10.5281/zenodo.18775726)
 
-<img src="man/figures/imgfile.png" width="200"> 
+<img src="man/figures/imgfile.png" alt="prolfqua hex sticker" width="200"> 
 
 # prolfqua - a comprehensive R package for Proteomics Differential Expression Analysis
 

_pkgdown.yml

@@ -14,3 +14,4 @@ navbar:
     github:
       icon: fa-github
       href: https://github.com/wolski/prolfqua
+      aria-label: GitHub repository