Contrast function mixes conditions in WaldTest output

[Roman-Si]

First of all, thank you for your very valuable package and thorough documentation

I encountered an unexpected behavior in the Contrasts function, specifically involving incorrect condition assignments during Wald test executions. In the sample data provided a protein is present in the conditions "Psfilm" and Psfoam" and absent in the condition "glucose" but the output of the Wald test shows this protein tested for the contrast "Psfoam - glucose".

To Reproduce I attach a sample file with 7 proteins, an annotation csv and an R script to reproduce.

Expected behavior The WaldTest output should write that the contrast is "Psfoam - Psfilm" for protein "g99273"

Screenshots

Desktop (please complete the following information):

• OS: Ubuntu 22.04.4 LTS
• prolfqua version: 1.2.3

Additional context Add any other context about the problem here. prolfqua-inputAnnotation.csv protein_intensities.csv prolfqua_contrast_testCode.txt

[wolski]

Thank you for sharing the reproducible example and spotting the problem. Yes it is unfortunately a bug in prolfqua.

Since you have 3 groups the model typically looks like this:

Call:
lm(formula = formula, data = x)

Coefficients:
     (Intercept)  condition_Psfilm  condition_Psfoam  
          23.743             1.847             2.131  

and the linear functions to compute the contasts, look like this:

> contr$get_linfct()
                     (Intercept) condition_Psfilm condition_Psfoam
Psfoam - Psfilm                0             -1.0              1.0
Psfoam - glucose               0              0.0              1.0
Psfilm - glucose               0              1.0              0.0

However for g99273 the modelling result is:

> modg99273
# A tibble: 1 × 9
# Groups:   protein_Id [1]
  protein_Id data             linear_model exists_lmer isSingular df.residual sigma nrcoef nrcoeff_not_NA
  <chr>      <list>           <list>       <lgl>       <lgl>            <dbl> <dbl>  <int>          <int>
1 g99273     <tibble [8 × 8]> <lm>         TRUE        TRUE                 1 0.361      2              2

and the linear model looks like this:

> modg99273$linear_model[1]
[[1]]

Call:
lm(formula = formula, data = x)

Coefficients:
     (Intercept)  condition_Psfoam  
          26.160            -1.043  

The intercept, in the model for g99273, estimates not the glucose group but the Psfilm group. Hence, we should have used a linear function, for Psfoam - Psfilm , i.e. c(0, 1) only.

An ad hoc fix would be:

mod <- prolfqua::build_model(
  lfqdata$data,
  formula_Condition,
  subject_Id = lfqdata$config$table$hierarchy_keys())

mod$modelDF <- mod$modelDF |> dplyr::filter(nrcoeff_not_NA == 3)

By this we remove all the models where not all 3 model parameter could not be estimated.

[wolski]

@Roman-Si I just pushed two commits that fix the error you have reported:

In the figure below, you can see that the contrast for protein g99273 is now correct.

Thank you for trying prolfqua and reporting problems.

Please update prolfqua

install.packages('remotes')
remotes::install_github('fgcz/prolfqua', dependencies = TRUE)

[wolski]

New release https://github.com/fgcz/prolfqua/releases/tag/v.1.2.4

[Roman-Si]

Hello, using another dataset with a 3x3 design I noticed again unexpected behaviour. The example protein here is present in 2 conditions but the groupAverage is used in all comparisons (even in alk_3d - alk_2d).

I attach sample code to reproduce and the initial peptide level data. I also noticed that TMP summarization creates artifacts where the intensities of some proteins are identical in all samples. This happens for proteins with one-sample peptides (like g004095, g012327 and others) since the row-median subtraction sets all rows to 0. This issue is also mentioned here when using median-polish for microarrays https://doi.org/10.1186/1471-2105-11-553. testcode_prolfqua.txt quantms_msstats_input.csv.gz

[wolski]

Thanks a lot for reporting the issue. You are observing two things. One is that we are getting group differences, although we do not have any observations in one of the groups. This is because you are using the ContrastMissing function, which estimates a LOD, substitutes the missing values with the LOD, and then computes the group averages. Then, you merged the contrast results using merge_contrasts_results. You now have two types of models in your data frame ( a linear model for the contrasts, which could be estimated using linear models, and the group average model for all those proteins with an excessive number of missing values). You can use the modelName column in the merged data.frame to see how each contrast was estimated.

mC <- ContrastsMissing$new(lfqdata = lfqdataNormalized, contrasts = Contrasts)
merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)$merged

The other problem is that we have standard error estimates that are 0, which is indeed caused by using Medpolish -> that the protein level estimates within the group are all equal (22.0). I do not like error estimates 0 (although running the contrasts through ContrastsModerated should shrink them to more reasonable estimates), and therefore , I updated the code for the contrast computation in ContrastsMissing, so that if std.error is 0, I am using the 75% Quantile of std.error of all proteins. I just pushed those changes.

You can use a different method for aggregation (sum_topN or lmrob). But those have their own issues.

Lastly, I suggest you apply the moderation to the merged contrasts:

merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)$merged
contr <- prolfqua::ContrastsModerated$new(merged)