In wolski/prolfqua: Proteomics Label Free Quantification

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
conflicted::conflict_prefer("filter", "dplyr")

Purpose

This vignette demonstrates how two conditions, e.g., treatment vs. control, can be compared and the differences statistically tested. We will again use the Ionstar dataset as an example of an LFQ experiment. This dataset was preprocessed with the MaxQuant software. After first examining the data using QC plots and then normalizing the data, we compare groups of replicates with the different dilutions. The output of the comparison is the difference in the mean intensities for quantified proteins (log2 fold-change) in each group, along with statistical parameters such as degrees of freedom, standard errors, p-value and the FDR.

Loading protein abundances from MaxQuant proteinGroups.txt

library(prolfqua)

Specify the path to the MaxQuant proteinGroups.txt file. The function tidyMQ_ProteinGroups will read the proteinGroups.txt file and convert it into a tidy table

xx <- prolfqua::sim_lfq_data_protein_config(Nprot = 100)
xx

Create the LFQData class instance and remove zeros from data (MaxQuant encodes missing values with zero).

lfqdata <- prolfqua::LFQData$new(xx$data, xx$config)
lfqdata$remove_small_intensities()
lfqdata$factors()

You can always convert the data into wide format.

lfqdata$to_wide()$data[1:3,1:7]

Visualization of not normalized data

After this first setting up of the analysis we show now how to normalize the proteins and the effect of normalization. Furthermore we use some functions to visualize the missing values in our data.

lfqplotter <- lfqdata$get_Plotter()
density_nn <- lfqplotter$intensity_distribution_density()

Visualization of missing data

lfqplotter$NA_heatmap()

lfqdata$get_Summariser()$plot_missingness_per_group()

lfqplotter$missigness_histogram()

Computing standard deviations, mean and CV.

Other important statistics such as coefficient of variation, means and standard deviations can be easily calculated using the get_Stats function and visualized with a violin plot.

stats <- lfqdata$get_Stats()
stats$violin()
prolfqua::table_facade( stats$stats_quantiles()$wide, paste0("quantile of ",stats$stat ))

stats$density_median()

Normalize protein intensities and show diagnostic plots

We normalize the data by $\log_2$ transforming and then $z-scaling$.

lt <- lfqdata$get_Transformer()
transformed <- lt$log2()$robscale()$lfq
transformed$config$table$is_response_transformed

pl <- transformed$get_Plotter()
density_norm <- pl$intensity_distribution_density()

gridExtra::grid.arrange(density_nn, density_norm)

pl$pairs_smooth()

p <- pl$heatmap_cor()

p

Fitting a linear model

For fitting linear models to the transformed intensities for all our proteins we have to first specify the model function and define the contrasts that we want to calculate.

transformed$factors()

formula_Condition <-  strategy_lm("transformedIntensity ~ group_")

# specify model definition
modelName  <- "Model"
Contrasts <- c("AvsC" = "group_A - group_Ctrl",
"BvsC" = "group_B - group_Ctrl")

Here we have to build the model for each protein.

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

In this plot we can see what factors in our model are mostly responsible for the adjusted p-values calculated from an analysis of variance.

mod$anova_histogram("FDR")

One also look what proteins do show different abundances in any of our five dilutions by looking at the FDR values of the analysis of variane.

aovtable <- mod$get_anova()
head(aovtable)
dim(aovtable)
xx <- aovtable |> dplyr::filter(FDR < 0.2)
signif <- transformed$get_copy()
signif$data <- signif$data |> dplyr::filter(protein_Id %in% xx$protein_Id)
hmSig <- signif$get_Plotter()$heatmap()

hmSig

Compute contrasts

Next we do calculate the statistics for our defined contrasts for all the proteins. For this we can use the Contrasts function.

contr <- prolfqua::Contrasts$new(mod, Contrasts)
v1 <- contr$get_Plotter()$volcano()

Alternatively, we can moderate the variance and using the Experimental Bayes method implemented in ContrastsModerated.

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

In the next figure it can be seen WEWinputNEEDED.

plotter <- contr$get_Plotter()
v2 <- plotter$volcano()
gridExtra::grid.arrange(v1$FDR,v2$FDR, ncol = 1)

plotter$ma_plotly()

#myProteinIDS <- c("sp|Q12246|LCB4_YEAST",  "sp|P38929|ATC2_YEAST",  "sp|Q99207|NOP14_YEAST")
myProteinIDS <- c("sp|P0AC33|FUMA_ECOLI",  "sp|P28635|METQ_ECOLI",  "sp|Q14C86|GAPD1_HUMAN")
dplyr::filter(contrdf, protein_Id %in% myProteinIDS)

Contrasts with missing value imputation

Use this method if there proteins with no observations in one of the groups. With the ContrastsMissing function, we can estimate difference in mean for proteins that are not observed in one group or condition. For this we are using the average expression at percentile 0.05 of the group where the protein is not quantified.

mC <- ContrastsMissing$new(lfqdata = transformed, contrasts = Contrasts)
colnames(mC$get_contrasts())

Finally we are merging the results and give priority to the results where we do not have missing values in one group.

merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)$merged
plotter <- merged$get_Plotter()
tmp <- plotter$volcano()
tmp$FDR

Look at proteins which could not be fitted using the linear model, if any.

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

moreProt <- transformed$get_copy()
moreProt$data <- moreProt$data |> dplyr::filter(protein_Id %in% merged$more$contrast_result$protein_Id)
moreProt$get_Plotter()$raster()

# here we do not get anything because there is nothing imputed!

GSEA Analyis

We can rank the proteins based on the log2FC or the t-statistic and subject them them to gene set enrichment analysis (see GSEA).

This example will run only if the following packages are installed on you machine:

• clusterProfiler (Bioconductor)
• org.Sc.sgd.db (Bioconductor)
• prora (github.com/protviz/prora)

#evalAll <- require("clusterProfiler") & require("org.Sc.sgd.db") & require("prora")
evalAll <- require("clusterProfiler") & require("org.Sc.sgd.db2") & require("prora")

library(clusterProfiler)
library(org.Sc.sgd.db)

bb <- prolfqua::get_UniprotID_from_fasta_header(merged$merged$get_contrasts(),
                                             idcolumn = "protein_Id")
bb <- prora::map_ids_uniprot(bb)
ranklist <- bb$statistic
names(ranklist) <- bb$P_ENTREZGENEID
res <- clusterProfiler::gseGO(
  sort(ranklist, decreasing = TRUE),
  OrgDb = org.Sc.sgd.db,
  ont = "ALL")

ridgeplot( res )

dotplot(res , showCategory = 30)

enrichplot::upsetplot(res)

The prolfqua package is described in [@Wolski2022.06.07.494524].

Session Info

sessionInfo()

References

wolski/prolfqua documentation built on Oct. 31, 2024, 9:22 a.m.