In wolski/prolfqua: Proteomics Label Free Quantification

knitr::opts_chunk$set(echo = TRUE)

alternatively use

dd <- prolfqua::sim_lfq_data_protein_config(Nprot = 100,weight_missing = 2)
dd$data$abundance |> is.na() |> sum()

Contrasts <- c("dilution.b-a" = "group_A - group_B", "dilution.c-e" = "group_A - group_Ctrl")
mh1 <- prolfqua::MissingHelpers$new(dd$data, dd$config, prob = 0.5,weighted = TRUE)
imputed <- mh1$get_contrasts(Contrasts = Contrasts)
mh2 <- prolfqua::MissingHelpers$new(dd$data, dd$config, prob = 0.5,weighted = FALSE)
imputed2 <- mh2$get_contrasts(Contrasts = Contrasts)
plot(imputed$estimate, imputed2$estimate)
abline(0 , 1 , col=2 , lwd=2)

mh1$get_LOD()
plot( imputed$estimate, -log10(imputed$p.value), pch = "*" )
points(imputed2$estimate, -log10(imputed2$p.value), col = 2, pch = "x")

modelling with linear models

Model with missing data

modelName <- "f_condtion_r_peptide"
formula_Protein <-
  prolfqua::strategy_lm("abundance  ~ group_",
              model_name = modelName)


mod <- prolfqua::build_model(
  dd$data,
  formula_Protein,
  modelName = modelName,
  subject_Id = dd$config$table$hierarchy_keys_depth())
mod$modelDF
mod$modelDF$nrcoeff_not_NA |> table()


mod$modelDF$isSingular |> table()
mod$modelDF |> nrow()
mod$get_anova()
prolfqua::model_summary(mod)

maxcoef <- max(mod$modelDF$nrcoeff_not_NA, na.rm = TRUE)
goodmods <- mod$modelDF |> dplyr::filter(isSingular == FALSE, exists_lmer == TRUE, nrcoeff_not_NA == maxcoef)

dim(goodmods)
xx <- lapply(goodmods$linear_model, vcov)
nr <- sapply(xx, nrow)
nr |> table()
nr <- sapply(xx, ncol)
nr |> table()
sum_matrix <- Reduce(`+`, xx)
sum_matrix/length(xx)

Model with lod imputation

loddata <- dd$data

loddata <- loddata |> dplyr::mutate(abundance = ifelse(is.na(abundance), mh1$get_LOD(), abundance))
modI <- prolfqua::build_model(
  loddata,
  formula_Protein,
  modelName = modelName,
  subject_Id = dd$config$table$hierarchy_keys_depth())


modI$modelDF$nrcoeff_not_NA |> table()
modI$modelDF$isSingular |> table()
modI$modelDF |> nrow()

allModels <- modI$modelDF$linear_model

xx <- lapply(allModels, vcov)
sum_matrix <- Reduce(`+`, xx)
sum_matrix/length(xx)

m <- (modI$modelDF$linear_model[[1]])
df.residual(m)
sigma(m)
vcov(m)

wolski/prolfqua documentation built on May 12, 2024, 10:16 p.m.