--- title: "Comparing Two Groups with prolfqua" author: "Witold E. Wolski" date: "`r Sys.Date()`" output: BiocStyle::html_document vignette: > %\VignetteIndexEntry{Comparing Two Groups with prolfqua} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console bibliography: - prolfqua.bib --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE) library(dplyr) ``` # Purpose This vignette demonstrates how two conditions, e.g., treatment vs. control, can be compared and the differences statistically tested. We use simulated protein abundances of an LFQ experiment with two treatment groups (`A`, `B`) and a control group (`Ctrl`). After first examining the data using QC plots and then normalizing the data, we compare each treatment group with the control. 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. # Simulated protein abundances ```{r createConfig} library(prolfqua) ``` `sim_lfq_data_protein_config()` simulates protein abundances in long format, one row per protein and sample, and returns them together with a matching `AnalysisConfiguration`. To read your own data, see the vignette [Creating configurations](https://fgcz.github.io/prolfqua/articles/CreatingConfigurations.html). ```{r LoadDataAndConfigure} xx <- prolfqua::sim_lfq_data_protein_config(Nprot = 100) xx ``` Create the `LFQData` class instance and remove zero or very small intensities, which most quantification software reports for missing values. ```{r removeSmallIntensities} lfqdata <- prolfqua::LFQData$new(xx$data, xx$config) lfqdata$remove_small_intensities() lfqdata$factors() ``` You can always convert the data into wide format. ```{r showWide, eval = TRUE} lfqdata$data_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. ```{r getPlotter} lfqplotter <- lfqdata$get_Plotter() density_nn <- lfqplotter$intensity_distribution_density() ``` ### Visualization of missing data ```{r makeMissingHeatmap, fig.cap="Heatmap of missing protein intensities, black - missing, white - present"} lfqplotter$na_heatmap() ``` ```{r missignessPerGroup, fig.cap="# of proteins with 0,1,...N missing values"} lfqdata$get_Summariser()$plot_missingness_per_group() ``` ```{r missignessHistogram, fig.cap="Intensity distribution of proteins depending on # of missing values"} 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. ```{r PlotCVDistributions, fig.cap="Violin plots of CVs in the different groups and among all groups"} stats <- lfqdata$get_Stats() stats$violin() prolfqua::table_facade( stats$stats_quantiles()$wide, paste0("quantile of ",stats$stat )) ``` ```{r plotCVsplitbyMedianIntensity, fig.cap="Distribution of CV's for top 50% and bottom 50% proteins by intensity."} stats$density_median() ``` ## Normalize protein intensities and show diagnostic plots We normalize the data by $\log_2$ transforming and then $z-scaling$. ```{r normalizedata} lt <- lfqdata$get_Transformer() transformed <- lt$log2()$robscale()$lfq transformed$is_transformed() ``` ```{r genplotNorm, fig.cap="Normalized intensities."} pl <- transformed$get_Plotter() density_norm <- pl$intensity_distribution_density() ``` ```{r showIntensityDistributions, fig.cap="Distribution of intensities before and after normalization."} gridExtra::grid.arrange(density_nn, density_norm) ``` ```{r plotScatterMatrix, fig.cap = "Scatterplot matrix"} pl$pairs_smooth() ``` ```{r createHeatmap} p <- pl$heatmap_cor() ``` ```{r plotHeatmap, fig.cap="Heatmap, Rows - proteins, Columns - samples", fig.align=5, fig.height=5} 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. ```{r lookatfactors} transformed$factors() ``` ```{r defineModelAndContrasts} formula_Condition <- strategy_lm("transformedIntensity ~ group_") # specify model definition modelName <- "Model" contr_spec <- c("AvsC" = "group_A - group_Ctrl", "BvsC" = "group_B - group_Ctrl") ``` Here we have to build the model for each protein. ```{r buildModel} mod <- prolfqua::build_model( transformed$data_long(), formula_Condition, subject_id = transformed$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. ```{r showANOVA, fig.cap="Distribution of adjusted p-values (FDR)"} mod$anova_histogram("FDR") ``` One can also look at which proteins show different abundances in any of the groups by looking at the FDR values of the analysis of variance. ```{r filterDataForFDR} aovtable <- mod$get_anova() head(aovtable) dim(aovtable) xx <- aovtable |> dplyr::filter(FDR < 0.2) signif <- transformed$get_copy() signif$set_data(signif$data_long() |> dplyr::filter(protein_Id %in% xx$protein_Id)) hmSig <- signif$get_Plotter()$heatmap() ``` ```{r showSigHeatmap, fig.cap="Heatmap for proteins with FDR < 0.2 in the analysis of variance"} 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. ```{r computeContrasts} contr <- prolfqua::Contrasts$new(mod, contr_spec) v1 <- contr$get_Plotter()$volcano() ``` Alternatively, we can moderate the variance using the empirical Bayes method implemented in `ContrastsModerated`. ```{r computeContrastsModerated} contr <- prolfqua::ContrastsModerated$new(contr) contrdf <- contr$get_contrasts() ``` Moderation shrinks the variance of each protein towards the variance of all proteins, which stabilises the p-values of proteins whose variance is estimated from few observations. The next figure compares the volcano plots without and with moderation. ```{r plotVolcanos, fig.cap="Volcano plot, Left panel - no moderation, Right panel - with moderation."} plotter <- contr$get_Plotter() v2 <- plotter$volcano() gridExtra::grid.arrange(v1$FDR,v2$FDR, ncol = 1) ``` ```{r showMAplot, fig.cap="MA plot showing the dependency of mean abundance with respect to the difference"} plotter$ma_plotly() ``` ```{r checkProteinsInOutput} myProteinIDS <- head(unique(contrdf$protein_Id), 3) dplyr::filter(contrdf, protein_Id %in% myProteinIDS) ``` # Contrasts with missing value imputation Use this method if there are proteins with no observations in one of the groups. The `lm_impute` facade refits the linear model for these proteins after imputing the missing observations with the limit of detection (LOD), and borrows the variance from the proteins that could be fitted. ```{r computeMissing, eval=TRUE} mC <- prolfqua::ContrastsLMImputeFacade$new(transformed, "~ group_", contr_spec)$contrast 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. ```{r mergeResults, fig.cap="Volcano plots for the two contrasts with missing value imputation from the lm_impute model."} 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. ```{r mergedMore} merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC) moreProt <- transformed$get_copy() moreProt$set_data(moreProt$data_long() |> 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! ``` # Gene Set Enrichment Analysis (GSEA) After computing contrasts, proteins can be ranked by log2 fold change or t-statistic and subjected to gene set enrichment analysis. The typical workflow is: 1. Extract UniProt IDs from the contrast results using `prolfqua::get_uniprot_id_from_fasta_header()` 2. Map UniProt IDs to Entrez Gene IDs using `AnnotationDbi::mapIds()` with a species-specific annotation package 3. Create a named ranked list (statistic values named by gene IDs) 4. Run GSEA using `clusterProfiler::gseGO()` or `clusterProfiler::gseKEGG()` For details see: - [clusterProfiler book](https://yulab-smu.top/biomedical-knowledge-mining-book/) -- comprehensive guide to enrichment analysis in R - [clusterProfiler Bioconductor page](https://bioconductor.org/packages/clusterProfiler/) - [AnnotationDbi](https://bioconductor.org/packages/AnnotationDbi/) -- for mapping between gene identifier types The `prolfqua` package is described in [@Wolski2022.06.07.494524]. # Session Info ```{r session_info} sessionInfo() ``` # References