Baldur Yeast Tutorial

In baldur: Bayesian Hierarchical Modeling for Label-Free Proteomics

1. Setup

This tutorial is quite fast and on a very simple data set (2 conditions only), for a more complicated tutorial on the setup please see vignette('baldur_ups_tutorial'). First we load baldur and setup the model dependent variables we need, then normalize the data and add the mean-variance trends.

library(baldur)
# Setup design matrix
yeast_design <- model.matrix(~0+factor(rep(1:2, each = 3)))
colnames(yeast_design) <- paste0('ng', c(50, 100))
# Compare the first and second column of the design matrix
# with the following contrast matrix
yeast_contrast <- matrix(c(-1, 1), nrow = 2)

# Set id column
id_col <- colnames(yeast)[1] # "identifier"

# Since baldur itself does not deal with missing data we remove the
# rows that have missing values for the purpose of the tutorial.
# Else, one would replace the filtering step with imputation but that is outside
# the scope of baldur
yeast_norm <- yeast %>%
  # Remove missing data
  tidyr::drop_na() %>%
  # Normalize data (this might already have been done if imputation was performed)
  psrn(id_col) %>%
  # Add mean-variance trends
  calculate_mean_sd_trends(yeast_design)

Importantly, note that the column names of the design matrix are unique subsets of the names of the columns within the conditions:

colnames(yeast)
#> [1] "identifier" "ng50_1"     "ng50_2"     "ng50_3"     "ng100_1"    "ng100_2"    "ng100_3"
colnames(yeast_design)
#> [1] "ng50"  "ng100"

This is essential for baldur to know which columns to use in calculations and to perform transformations.

2. Mean-Variance trends and Gamma Regression fitting

Next is to infer the mixture of the data and to estimate the parameters needed for baldur. First we will setup the needed variables for using baldur without partitioning the data. Then, partitioning and setting up baldur after trend-partitioning

# Fit the gamma regression
gr_model <- fit_gamma_regression(yeast_norm, sd ~ mean)
# Estimate the uncertainty
unc_gr <- estimate_uncertainty(gr_model, yeast_norm, id_col, yeast_design)

3. Run the sampling procedure

Finally we sample the posterior of each row in the data. First we sample assuming a single trend, then using the partitioning.

# Single trend
gr_results <- gr_model %>%
  # Add hyper-priors for sigma
  estimate_gamma_hyperparameters(yeast_norm) %>%
  infer_data_and_decision_model(
    id_col,
    yeast_design,
    yeast_contrast,
    unc_gr,
    clusters = 10 # I highly recommend using parallel workers/clusters
  )               # this will greatly reduce the speed of running baldur
# The top hits then looks as follows:
gr_results %>%
  dplyr::arrange(err)
#> # A tibble: 1,802 × 22
#>    identifier         comparison       err   lfc lfc_025 lfc_50 lfc_975 lfc_eff lfc_rhat  sigma sigma_025 sigma_50 sigma_975 sigma_eff
#>    <chr>              <chr>          <dbl> <dbl>   <dbl>  <dbl>   <dbl>   <dbl>    <dbl>  <dbl>     <dbl>    <dbl>     <dbl>     <dbl>
#>  1 Cre09.g406050.t1.… ng100 vs … 7.63e-224  6.18   5.79    6.18    6.59   2079.     1.00 0.124     0.0676   0.114     0.243      1372.
#>  2 sp|P37302|APE3_YE… ng100 vs … 3.08e-184  1.51   1.41    1.51    1.61   3731.     1.00 0.0450    0.0245   0.0415    0.0869     1403.
#>  3 Cre12.g554100.t1.… ng100 vs … 1.88e-176  1.61   1.49    1.61    1.72   3203.     1.00 0.0493    0.0265   0.0451    0.0987     1072.
#>  4 sp|P38788|SSZ1_YE… ng100 vs … 5.72e-151  1.07   0.991   1.07    1.15   3696.     1.00 0.0357    0.0192   0.0324    0.0721     1154.
#>  5 Cre14.g616700.t1.… ng100 vs … 4.83e-118 -4.55  -4.96   -4.55   -4.14   2574.     1.00 0.148     0.0799   0.134     0.297       893.
#>  6 Cre10.g420750.t1.… ng100 vs … 7.84e-113  4.16   3.80    4.16    4.52   3056.     1.00 0.153     0.0847   0.140     0.305      1290.
#>  7 Cre06.g306550.t1.… ng100 vs … 4.38e-102  4.20   3.82    4.20    4.60   2972.     1.00 0.148     0.0816   0.137     0.287      1445.
#>  8 Cre12.g533100.t1.… ng100 vs … 9.87e- 98  1.41   1.29    1.41    1.55   3087.     1.00 0.0593    0.0325   0.0539    0.119      1232.
#>  9 sp|P09624|DLDH_YE… ng100 vs … 1.49e- 97  1.48   1.34    1.48    1.62   3135.     1.00 0.0592    0.0326   0.0539    0.119      1422.
#> 10 sp|P32481|IF2G_YE… ng100 vs … 1.38e- 90  1.59   1.43    1.59    1.74   3145.     1.00 0.0680    0.0370   0.0624    0.135      1500.
#> # ℹ 1,792 more rows
#> # ℹ 8 more variables: sigma_rhat <dbl>, lp <dbl>, lp_025 <dbl>, lp_50 <dbl>, lp_975 <dbl>, lp_eff <dbl>, lp_rhat <dbl>,
#> #   warnings <list>

Here err is the probability of error, i.e., the two tail-density supporting the null-hypothesis, lfc is the estimated log$_2$-fold change, sigma is the common variance, and lp is the log-posterior. Columns without suffix shows the mean estimate from the posterior, while the suffixes _025, _50, and _975, are the 2.5, 50.0, and 97.5, percentiles, respectively. The suffixes _eff and _rhat are the diagnostic variables returned by rstan (please see the Stan manual for details). In general, a larger _eff indicates a better sampling efficiency, and _rhat compares the mixing within chains against between the chains and should be smaller than 1.05.

4. Running Baldur with Latent Gamma Mixture Regression estimation

First we fit the LGMR model:

yeast_lgmr <- fit_lgmr(yeast_norm, id_col, lgmr_model, cores = 5)

We can print the model with print and extract parameters of interest with coef:

print(yeast_lgmr, pars = c("coef", "aux"))
#> 
#> LGMR Model
#>  μ = exp(-1.848 - 0.326 f(ȳ)) + κ exp(θ(7.521 - 0.474 f(ȳ)))
#> 
#>  Auxiliary:
#>       mean   se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff  Rhat
#> α     4.48  0.007271  0.299  3.912  4.283  4.477  4.680  5.093   1689     1
#> RMSE  0.52  0.000389  0.025  0.475  0.503  0.518  0.536  0.573   4123     1
#> 
#> 
#>  Coefficients:
#>         mean   se_mean      sd    2.5%     25%     50%     75%   97.5%  n_eff  Rhat
#> γ_0    7.521  0.000729  0.0505   7.424   7.487   7.521   7.555   7.622   4795     1
#> γ_0L  -1.848  0.000699  0.0299  -1.905  -1.868  -1.848  -1.828  -1.788   1828     1
#> γ_ȳ    0.326  0.000362  0.0246   0.279   0.309   0.325   0.342   0.374   4606     1
#> γ_ȳL   0.474  0.000546  0.0442   0.387   0.444   0.474   0.504   0.561   6562     1
# Extract the regression, alpha, and theta parameters and the NRMSE.
yeast_lgmr_coef <- coef(yeast_lgmr, pars = "all")

Baldur allows for two ways to plot the LGMR model, plot_lgmr_regression, and plot_regression_field. The first plots lines of three cases of $\theta$, 0, 0.5, and 1, and colors each peptide according to their infered $\theta$. They can be plotted accordingly:

plot_lgmr_regression(yeast_lgmr)
plot_regression_field(yeast_lgmr, rng = 25)

We can then estimate the uncertainty similar to the GR case:

unc_lgmr <- estimate_uncertainty(yeast_lgmr, yeast_norm, id_col, yeast_design)

Then running the data and decision model:

# Single trend
lgmr_results <- yeast_lgmr %>%
  # Add hyper-priors for sigma
  estimate_gamma_hyperparameters(yeast_norm, id_col) %>%
  infer_data_and_decision_model(
    id_col,
    yeast_design,
    yeast_contrast,
    unc_lgmr,
    clusters = 10
  )

5. Visualization of the results

baldur have two ways of visualizing the results 1) plotting sigma vs LFC and 2) Volcano plots. To plot sigma against LFC we use plot_sa:

gr_results %>%
  plot_sa(
    alpha = .05, # Level of significance
    lfc = 1      # Add LFC lines
  )

lgmr_results %>%
  plot_sa(
    alpha = .05, # Level of significance
    lfc = 1      # Add LFC lines
  )

While it is hard to see with this few examples, in general a good decision is indicated by a lack of a trend between $\sigma$ and LFC. To make a volcano plot one uses plot_volcano in a similar fashion to plot_sa:

gr_results %>%
  plot_volcano(
    alpha = .05, # Level of significance
    lfc = 1      # Add LFC lines
  )

lgmr_results %>%
  plot_volcano(
    alpha = .05, # Level of significance
    lfc = 1      # Add LFC lines
  )

baldur documentation built on Sept. 18, 2023, 9:07 a.m.