Skip to contents

Build a contrast analysis using one of several statistical methods

Source: R/build_contrast_analysis.R
build_contrast_analysis.Rd

The factory of prolfqua's modelling: method selects a facade class from the facade registry, and the facade runs the full pipeline from strategy construction through modelling to contrast computation. Every facade returns its results through the same interface (get_contrasts(), get_Plotter(), get_missing(), to_wide()), so the results of different methods have the same columns and can be compared or combined.

Usage

build_contrast_analysis(lfqdata, modelstr, contrasts, method = "lm", ...)

Arguments

lfqdata

an LFQData object

modelstr

model formula string without the response variable (e.g. "~ group_"). The response is taken automatically from lfqdata$get_config()$get_response().

contrasts

named character vector of contrasts (e.g. c("A_vs_B" = "group_A - group_B"))

method

a registered facade key. The built-in keys are: linear models "lm", "rlm", "rfit"; limma "limma", "limma_voom", "deqms", "deqms_voom", "limpa", "limpa_nested"; limit-of-detection imputation "lm_impute", "limma_impute", "limma_voom_impute", "rfit_impute"; peptide-level models "lmer_nested", "ropeca_nested"; detection models "firth", "firth_nested", "binomial_nested", "binomial_firth_nested"; the hurdle model "hurdle_nested", combining intensity and detection. Downstream packages may add more via register_facade. The authoritative list is names(list_facades()). Defaults to "lm".

...

additional arguments forwarded to the underlying strategy function (e.g. trend, robust for strategy_limma)

Value

a facade object of the class registered for method, see lookup_facade.

Details

Methods ending in _nested take peptide- or precursor-level input and report protein-level contrasts; the others take input at the level they report, usually aggregated proteins. "limpa" needs input aggregated with AggregateLimpa.

See also

Other modelling: AnovaExtractor, Contrasts, ContrastsDEqMSFacade, ContrastsDEqMSVoomFacade, ContrastsFacadeBase, ContrastsFirth, ContrastsFirthFacade, ContrastsFirthNestedFacade, ContrastsHurdleFacade, ContrastsHurdleNestedFacade, ContrastsLMFacade, ContrastsLMImputeFacade, ContrastsLimma, ContrastsLimmaFacade, ContrastsLimmaImputeFacade, ContrastsLimmaVoomFacade, ContrastsLimmaVoomImputeFacade, ContrastsLimpaFacade, ContrastsLimpaNestedFacade, ContrastsLmerNestedFacade, ContrastsModerated, ContrastsModeratedDEqMS, ContrastsPlotter, ContrastsRLMFacade, ContrastsROPECA, ContrastsROPECANestedFacade, ContrastsRfitFacade, ContrastsRfitImputeFacade, ContrastsTable, INTERNAL_FUNCTIONS_BY_FAMILY, LR_test(), Model, ModelFirth, ModelLimma, StrategyLM, StrategyLimma, StrategyLimpa, StrategyLmer, StrategyLogistf, StrategyRLM, StrategyRfit, build_model(), build_model_glm_peptide(), build_model_glm_protein(), build_model_impute(), build_model_limma(), build_model_limma_impute(), build_model_limma_voom(), build_model_limma_voom_impute(), build_model_limpa(), build_model_logistf(), compute_borrowed_variance(), compute_borrowed_variance_limma(), compute_contrast(), compute_lmer_contrast(), df.residual.rfit_prolfqua(), get_anova_df(), get_complete_model_fit(), get_p_values_pbeta(), group_label(), impute_from_model(), impute_refit_singular(), linfct_from_model(), linfct_matrix_contrasts(), list_facades(), lookup_facade(), merge_contrasts_results(), model_analyse(), model_summary(), moderated_p_deqms(), moderated_p_deqms_long(), moderated_p_limma(), new_imputed_model(), pivot_model_contrasts_to_wide(), register_facade(), sigma.rfit_prolfqua(), sim_build_models_lm(), sim_build_models_lmer(), sim_build_models_logistf(), sim_make_model_lm(), sim_make_model_lmer(), strategy_limma(), strategy_limpa(), strategy_logistf(), summary_ROPECA_median_p_scaled(), vcov.rfit_prolfqua()

Examples

istar <- sim_lfq_data_protein_config(Nprot = 20)
#> creating sampleName from file_name column
#> completing cases
#> completing cases done
#> setup done
lfqdata <- LFQData$new(istar$data, istar$config)
lfqdata$rename_response("transformedIntensity")
contrasts <- c("A_vs_Ctrl" = "group_A - group_Ctrl")

fa_lm    <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "lm")
head(fa_lm$get_contrasts())
#> determine linear functions:
#> get_contrasts -> contrasts_linfct
#> contrasts_linfct
#> Joining with `by = join_by(protein_Id, contrast)`
#> # A tibble: 6 × 16
#>   modelName estimate_type protein_Id  contrast     diff avgAbd
#>   <chr>     <chr>         <chr>       <chr>       <dbl>  <dbl>
#> 1 lm        observed      0EfVhX~5954 A_vs_Ctrl  2.72     23.2
#> 2 lm        observed      0m5WN4~1448 A_vs_Ctrl  0.600    17.4
#> 3 lm        observed      7cbcrd~8305 A_vs_Ctrl  2.59     27.0
#> 4 lm        observed      9VUkAq~4562 A_vs_Ctrl  0.0679   19.4
#> 5 lm        observed      At886V~3296 A_vs_Ctrl -1.01     19.1
#> 6 lm        observed      BEJI92~9143 A_vs_Ctrl -0.873    20.9
#> # ℹ 10 more variables: std.error.unmoderated <dbl>, df.unmoderated <int>,
#> #   std.error <dbl>, statistic <dbl>, df <dbl>, p.value <dbl>, conf.low <dbl>,
#> #   conf.high <dbl>, sigma <dbl>, FDR <dbl>

fa_limma <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "limma")
head(fa_limma$get_contrasts())
#> # A tibble: 6 × 16
#>   modelName estimate_type protein_Id  contrast     diff    FDR
#>   <chr>     <chr>         <chr>       <chr>       <dbl>  <dbl>
#> 1 limma     observed      0EfVhX~5954 A_vs_Ctrl  2.72   0.188 
#> 2 limma     observed      0m5WN4~1448 A_vs_Ctrl  0.600  0.623 
#> 3 limma     observed      7cbcrd~8305 A_vs_Ctrl  2.59   0.0271
#> 4 limma     observed      9VUkAq~4562 A_vs_Ctrl  0.0679 0.967 
#> 5 limma     observed      At886V~3296 A_vs_Ctrl -1.01   0.623 
#> 6 limma     observed      BEJI92~9143 A_vs_Ctrl -0.873  0.623 
#> # ℹ 10 more variables: std.error.unmoderated <dbl>, df.unmoderated <dbl>,
#> #   std.error <dbl>, statistic <dbl>, p.value <dbl>, sigma <dbl>, df <dbl>,
#> #   conf.low <dbl>, conf.high <dbl>, avgAbd <dbl>

fa_deqms <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "deqms")
head(fa_deqms$get_contrasts())
#> determine linear functions:
#> get_contrasts -> contrasts_linfct
#> contrasts_linfct
#> Joining with `by = join_by(protein_Id, contrast)`
#> # A tibble: 6 × 16
#>   modelName estimate_type contrast  protein_Id     diff avgAbd
#>   <chr>     <chr>         <chr>     <chr>         <dbl>  <dbl>
#> 1 deqms     observed      A_vs_Ctrl 0EfVhX~5954  2.72     23.2
#> 2 deqms     observed      A_vs_Ctrl 0m5WN4~1448  0.600    17.4
#> 3 deqms     observed      A_vs_Ctrl 7cbcrd~8305  2.59     27.0
#> 4 deqms     observed      A_vs_Ctrl 9VUkAq~4562  0.0679   19.4
#> 5 deqms     observed      A_vs_Ctrl At886V~3296 -1.01     19.1
#> 6 deqms     observed      A_vs_Ctrl BEJI92~9143 -0.873    20.9
#> # ℹ 10 more variables: std.error.unmoderated <dbl>, df.unmoderated <int>,
#> #   std.error <dbl>, statistic <dbl>, df <int>, p.value <dbl>, conf.low <dbl>,
#> #   conf.high <dbl>, sigma <dbl>, FDR <dbl>

istar_pep <- sim_lfq_data_peptide_config()
#> creating sampleName from file_name column
#> completing cases
#> completing cases done
#> setup done
lfqdata_pep <- LFQData$new(istar_pep$data, istar_pep$config)
lfqdata_pep <- lfqdata_pep$get_Transformer()$log2()$lfq
#> Column added : log2_abundance

fa_lmer <- build_contrast_analysis(
  lfqdata_pep,
  "~ group_ + (1 | peptide_Id) + (1 | sampleName)",
  contrasts,
  method = "lmer_nested"
)
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning: There were 4 warnings in `dplyr::mutate()`.
#> The first warning was:
#> ℹ In argument: `linear_model = purrr::map(data, model_strategy$model_fun, pb =
#>   pb)`.
#> ℹ In group 2: `protein_Id = "7cbcrd~5725"`.
#> Caused by warning:
#> ! grouping factors must have > 1 sampled level
#> ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
head(fa_lmer$get_contrasts())
#> determine linear functions:
#> get_contrasts -> contrasts_linfct
#> contrasts_linfct
#> Joining with `by = join_by(protein_Id, contrast)`
#> # A tibble: 6 × 16
#>   modelName   estimate_type protein_Id  contrast       diff avgAbd
#>   <chr>       <chr>         <chr>       <chr>         <dbl>  <dbl>
#> 1 lmer_nested observed      0EfVhX~0087 A_vs_Ctrl -0.000832   4.34
#> 2 lmer_nested observed      BEJI92~5282 A_vs_Ctrl  0.322      4.22
#> 3 lmer_nested observed      Fl4JiV~8625 A_vs_Ctrl -0.0413     4.38
#> 4 lmer_nested observed      HvIpHG~9079 A_vs_Ctrl -0.372      4.40
#> 5 lmer_nested observed      JcKVfU~9653 A_vs_Ctrl -0.107      5.05
#> 6 lmer_nested observed      SGIVBl~5782 A_vs_Ctrl  0.0307     4.68
#> # ℹ 10 more variables: std.error.unmoderated <dbl>, df.unmoderated <dbl>,
#> #   std.error <dbl>, statistic <dbl>, df <dbl>, p.value <dbl>, conf.low <dbl>,
#> #   conf.high <dbl>, sigma <dbl>, FDR <dbl>

fa_ropeca <- build_contrast_analysis(lfqdata_pep, "~ group_", contrasts, method = "ropeca_nested")
head(fa_ropeca$get_contrasts())
#> determine linear functions:
#> get_contrasts -> contrasts_linfct
#> contrasts_linfct
#> Joining with `by = join_by(protein_Id, peptide_Id, contrast)`
#> # A tibble: 6 × 14
#> # Groups:   contrast [1]
#>   modelName   estimate_type protein_Id contrast avgAbd    diff     FDR statistic
#>   <chr>       <chr>         <chr>      <chr>     <dbl>   <dbl>   <dbl>     <dbl>
#> 1 ropeca_nes… observed      0EfVhX~00… A_vs_Ct…   4.27 -0.0742 5.28e-2     -1.75
#> 2 ropeca_nes… observed      7cbcrd~57… A_vs_Ct…   4.51  0.741  9.91e-5      8.79
#> 3 ropeca_nes… observed      9VUkAq~47… A_vs_Ct…   4.47 -0.598  6.91e-6    -12.7 
#> 4 ropeca_nes… observed      BEJI92~52… A_vs_Ct…   4.23  0.277  1.87e-3      3.94
#> 5 ropeca_nes… observed      CGzoYe~21… A_vs_Ct…   4.76 -0.310  3.74e-5     -9.26
#> 6 ropeca_nes… observed      DoWup2~58… A_vs_Ct…   4.43  0.295  1.38e-6     14.7 
#> # ℹ 6 more variables: std.error <dbl>, df <int>, p.value <dbl>, conf.low <dbl>,
#> #   conf.high <dbl>, sigma <dbl>

fa_firth <- build_contrast_analysis(lfqdata, "~ group_", contrasts, method = "firth")
#> completing cases
#> Joining with `by = join_by(protein_Id)`
head(fa_firth$get_contrasts())
#> determine linear functions:
#> get_contrasts -> contrasts_linfct
#> contrasts_linfct
#> Joining with `by = join_by(protein_Id, contrast)`
#> # A tibble: 6 × 16
#> # Groups:   contrast [1]
#>   modelName estimate_type protein_Id  contrast  sigma    df      diff   FDR
#>   <chr>     <chr>         <chr>       <chr>     <dbl> <int>     <dbl> <dbl>
#> 1 firth     observed      0EfVhX~5954 A_vs_Ctrl     1     9  1.07e-15     1
#> 2 firth     observed      0m5WN4~1448 A_vs_Ctrl     1     9 -2.20e+ 0     1
#> 3 firth     observed      7cbcrd~8305 A_vs_Ctrl     1     9  1.07e-15     1
#> 4 firth     observed      9VUkAq~4562 A_vs_Ctrl     1     9 -1.35e+ 0     1
#> 5 firth     observed      At886V~3296 A_vs_Ctrl     1     9  5.58e-16     1
#> 6 firth     observed      BEJI92~9143 A_vs_Ctrl     1     9 -1.35e+ 0     1
#> # ℹ 8 more variables: std.error <dbl>, statistic <dbl>, p.value <dbl>,
#> #   conf.low <dbl>, conf.high <dbl>, avgAbd <dbl>, std.error.unmoderated <dbl>,
#> #   df.unmoderated <int>