In wolski/prolfqua: Proteomics Label Free Quantification

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

Purpose

In this tutorial, we delve into the concept of using multiple factors, also known as explanatory variables, to model the observed variance in your data. We will demonstrate this by modeling data with two factors and their interaction.

Examples of data where two explanatory variables are needed to explain the variance in the data are for instance: - Two cell lines (X) and (Z), for each of which we measured a control condition (A) and a treatment condition (B). - An experiment where samples from a control condition (A) and treatment condition (B) were measured in two batches, X and Y, and there is a batch effect we must account for. - A combination of treatments A and B results in factors such as FA with levels placeboA and A and FB with levels placeboB and B.

Let's assume that the underlying dataset is generated in a course held annually. The context is that yeast is grown on glucose in one condition (A), and in the other condition (B), yeast is grown on glycerol and ethanol. Here, we are looking into the results of two different batches (X and Z), where other people performed the wet lab work, and even different LC-MS instruments were involved. It is, therefore, essential to model the batch variable for these two similar datasets.

We are also modeling the interaction between the two explanatory variables batch and condition for demonstration purposes. In this case, having a significant interaction term would mean the protein is expressed more in the Glucose condition in one batch. In contrast, the same protein is more abundant in the Ethanol condition in the other batch.

An in depth introduction to modelling and testing interactions using linear models can be found here.

Model Fitting

We use simulated data generated using the function sim_lfq_data_2Factor_config. Interesting here is the definition of the model. If interaction shall be included in the model a asterix should be used while if no interaction should be taken into account a plus should be used in the model definition. Also we can directly specify what comparisons we are interested in by specifying the respective contrasts.

conflicted::conflict_prefer("filter", "dplyr")

#data_Yeast2Factor <- prolfqua::prolfqua_data("data_Yeast2Factor")
data_2Factor <- prolfqua::sim_lfq_data_2Factor_config(
  Nprot = 200,
  with_missing = TRUE,
  weight_missing = 2)
data_2Factor <- prolfqua::LFQData$new(data_2Factor$data, data_2Factor$config)

pMerged <- data_2Factor
pMerged$factors()

formula_Batches <-
  prolfqua::strategy_lm("abundance ~ Treatment * Background ")

# specify model definition
Contrasts <- c("TA - TB" = "TreatmentA - TreatmentB",
               "BX - BY" = "BackgroundX - BackgroundZ",
               "AvsB_gv_BackgroundX" = "`TreatmentA:BackgroundX` - `TreatmentB:BackgroundX`",
               "AvsB_gv_BackgroundZ" = "`TreatmentA:BackgroundZ` - `TreatmentB:BackgroundZ`",
               "Interaction" = "AvsB_gv_BackgroundX - AvsB_gv_BackgroundZ")

We are then building our model as we specified it before for each protein.

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

Now, we can visualize the effect of our factors that are specified in the model. This indicates in an elegant way what factors are the most important ones.

mod$anova_histogram()$plot

ANOVA

To examine proteins with a significant interaction between the two factors treatment and batch filtering for the factor condition_:batch_.

ANOVA <- mod$get_anova()
ANOVA |> dplyr::filter(factor == "Treatment:Background") |> dplyr::arrange(FDR) |> head(5)
protIntSig <- ANOVA |> dplyr::filter(factor == "Treatment:Background") |>
  dplyr::filter(FDR < 0.01)
protInt <-  pMerged$get_copy()
protInt$data <- protInt$data[protInt$data$protein_Id %in% protIntSig$protein_Id[1:6],]

These proteins can easily be visualized using the boxplot function from the plotter objects in prolfqua

ggpubr::ggarrange(plotlist = protInt$get_Plotter()$boxplots()$boxplot)

Compute and analyse with the specified contrasts

Next, we want to calculate the statistical results for our group comparisons that have been specified in our contrast definition. Here we are using the moderated statistics which implements the concept of pooled variance for all proteins.

contr <- prolfqua::ContrastsModerated$new(prolfqua::Contrasts$new(mod, Contrasts))
contrdf <- contr$get_contrasts()

These results can be visualized with e.g a volcano or a MA plot.

plotter <- contr$get_Plotter()
plotter$volcano()$FDR
plotter$ma_plot()

Analyse contrasts with missing data imputation

Still using the approach above, we can only estimate group averages in case there is at least one measurement for each protein in each group/condition. In the case of missing data for one condition, we can use the ContrastsMissing function where the 10th percentile expression of all proteins is used for the estimate of the missing condition.

contrSimple <- prolfqua::ContrastsMissing$new(pMerged, Contrasts)
contrdfSimple <- contrSimple$get_contrasts()
pl <- contrSimple$get_Plotter()
pl$histogram_diff()
pl$volcano()$FDR

Merge nonimputed and imputed data.

For the moderated model, we can only get results if we have enough valid data points. With the group average model we can get group estimates for all proteins. Therefore, we are merging the statistics for both approaches while we are preferring the values of the moderated model. Also these results can again be visualized in a volcano plot.

dim(contr$get_contrasts())
dim(contrSimple$get_contrasts())

mergedContrasts <- prolfqua::merge_contrasts_results(prefer = contr, add = contrSimple)$merged
cM <- mergedContrasts$get_Plotter()
plot <- cM$volcano()
plot$FDR

Look at Proteins with significant interaction term.

sigInteraction <- mergedContrasts$contrast_result |> 
  dplyr::filter(contrast == "Interaction" & FDR < 0.2)

protInt <-  pMerged$get_copy()
protInt$data <- protInt$data[protInt$data$protein_Id %in% sigInteraction$protein_Id,]

protInt$get_Plotter()$raster()

hm <- protInt$get_Plotter()$heatmap()

hm

Alternative model specification

We compute the same contrasts as above but using only one factor and subgroups "A_X", "A_Z", "B_X", "B_Z".

We start by simulating the data.

data_1Factor <- prolfqua::sim_lfq_data_2Factor_config(
  Nprot = 200,
  with_missing = TRUE,
  weight_missing = 2, TWO = FALSE)
data_1Factor <- prolfqua::LFQData$new(data_1Factor$data, data_1Factor$config)


data_1Factor$response()

Instead of two factors we now have one factor Group with four levels r data_1Factor$factors()$Group |> table().

knitr::kable(data_1Factor$factors())

We specify the model formula and the same contrasts as for the two factor model but using only one factor and the subgroups.

formula_Batches <-
  prolfqua::strategy_lm("abundance ~ Group")

# specify model definition
Contrasts <- c("TA - TB" = "(GroupA_X + GroupA_Z)/2 - (GroupB_X + GroupB_Z)/2",
               "BX - BY" = "(GroupA_X + GroupB_X)/2 - (GroupA_Z + GroupB_Z)/2",
               "AvsB_gv_BackgroundX" = "GroupA_X - GroupB_X",
               "AvsB_gv_BackgroundZ" = "GroupA_Z - GroupB_Z",
               "Interaction" = "AvsB_gv_BackgroundX - AvsB_gv_BackgroundZ")

mod <- prolfqua::build_model(data_1Factor$data, formula_Batches,
  subject_Id = pMerged$config$table$hierarchy_keys() )

contr <- prolfqua::ContrastsModerated$new(prolfqua::Contrasts$new(mod, Contrasts))
contrdfONE <- contr$get_contrasts()

We now compare the contrasts computed from the model with two factors with those obtained from the model with one factor. We can see that the contrast estimates for difference, t-statistics, p.value and FDR are the same.

xx <- dplyr::inner_join(contrdf , contrdfONE, by = c("protein_Id","contrast"), suffix = c(".TWO",".ONE"))
par(mfrow = c(2,2))
plot(xx$diff.ONE, xx$diff.TWO)
plot(xx$statistic.ONE, xx$statistic.TWO)
plot(xx$FDR.ONE, xx$FDR.TWO)
plot(xx$p.value.ONE, xx$p.value.TWO)

Likelihood ratio Test for models with more factors

In cases where you have more then one factor possibly explaining the variance in your data, you can use the likelihood ratio test, to examine which factor to include into the statistical model. For more details see the LR_test function documentation and example code. (To open the documentation run ?LR_test in the R console.)

Testing interaction computation

url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/spider_wolff_gorb_2013.csv"
filename <- "spider_wolff_gorb_2013.csv"
downloader::download(url, filename)
spider <- read.csv(filename, skip = 1)
XA <- lm(friction ~ type * leg, data = spider)
summary(XA)

head(spider)
spider <- spider |> tidyr::unite(legtype , leg, type, remove = FALSE)
XF <- lm(friction ~ legtype - 1, data = spider)
library(multcomp)
summary(glht(XF, linfct = rbind(i1 = c(1, -1, -1, 1 ,0, 0, 0, 0))))
summary(glht(XF, linfct = rbind(i1 = c(1, -1, 0, 0, -1, 1, 0, 0 ))))
summary(glht(XF, linfct = rbind(i1 = c(1, -1, 0, 0,  0, 0, -1, 1  ))))

Session Info

sessionInfo()

References

wolski/prolfqua documentation built on Sept. 7, 2024, 3:06 a.m.