Biostatistics - Basics of experimental design, data processing and differential analysis

Author

Affiliation

Lieven Clement

Ghent University

1 Motivating example (Brunner et al. 2022)

Brunner et al. 2022. Mol Syst Biol. 18(3): e10798. doi: 10.15252/msb.202110798

• Single cell study on 231 HeLa cells with drug induced cell cycle arrest for which over 2500 proteins were measured with MS-SCP.

• Authors report:

“The proteomes of the different cell cycle states grouped together in a principal component analysis (PCA) plot”.

“Our single‐cell data set also highlighted proteins not previously associated with the cell cycle and the G2/M transition”

1.1 Reanalysis of the results

We could reproduce the results:

Another way to look at the same plot

What happened

Confounding between acquisition batch and cell cycle arrest!!

To consult the statistician after an experiment is finished is often merely to ask him to conduct a post mortem examination. He can perhaps say what the experiment died of. – Ronald Fisher

Thinking about design before your experiment is key!!!!

1.2 Why stats

• Statistics is the science on
  1. collecting (experimental design),
  2. exploring (data exploration) and
  3. learning from data and to generalize what we observe in the sample towards the population while quantifying, controlling and reporting variability and uncertainty (statistical modelling and statistical inference).

1.3 Structure of the lecture

1. Understand properties of MS-based proteomics data / Data Processing
2. Differential Abundance analysis of MS-based proteomics Data / Statistical Modeling & Inference
3. Some important aspects of experimental design

2 Data Processing

2.1 Data Background

2.1.1 Data acquisition

Data Dependent Acquisition

Overview of an LFQ-based proteomics workflow.

• Peptide Characteristics

  • Modifications
  • Ionisation Efficiency: huge variability
  • Identification
    • Misidentification \(\rightarrow\) outliers
    • MS\(^2\) selection on peptide abundance
    • Context depending missingness
    • Non-random missingness
    • Random missingness

\(\rightarrow\) Unbalanced pepide identifications across samples and messy data

Data Independent Acquisition

Data-independent acquisition.

• m/z range into sequential windows (e.g., 400–425, 425–450 m/z)

• MS2 spectra are complex as many peptides mixed

• Deconvolution of the MS2 signal

• With dedicated software such as Spectronaut or DIA-NN

  • Identification based on spectral libraries
  • Library free: FASTA database to predict in silico spectra, retention times, and ion mobilities, which are then used to search the DIA data
  • Use of decoys sequences to estimate false discovery rate of ID

• Quantification using MS2 and/or MS1 peaks

2.1.2 Level of quantification

• MS-based proteomics returns peptides or precursors: pieces of proteins

• Quantification commonly required on the protein level

2.1.3 Experimental context

• In-house DIA spike-in study to illustrate the analysis (subset of (Staes et al. 2024)).

• spiked digested UPS proteins in yeast background at the following ratio’s \(\rightarrow\) yeast:ups ratio 10:1, 10:2, 10:4, 10:8, 10:10

• Each sample was analyzed in triplicate

• Here use the data from the middle 3 spike-in ratio’s (10:2, 10:4 and 10:8)
• Searched using DIA-NN 2.2.0.

2.2 Precursor Filtering

• Remove questionable identifications: Low ID Q-value, Decoys, Contaminants, …
• Only keep Proteotypic precursors
• Remove precursors/PSMs with many missing values
• …

Bourgon et al. (2010) showed that filtering is only appropriate if the statistic used for filtering is independent of the statistic used for hypothesis testing under the null hypothesis.

2.3 Log-transformation

• Precursor: ADSRDPASDQMQHWK, charge-state 3 from UPS|P04040ups|CATA_HUMAN_UPS

Plot upon log2 -transformation

• Mean-variance relation \(\rightarrow\) Solved by log2 tranformation
• Other property of log transformation:

\[log_2FC_{b-a} = log_2b - log_2a = log_2 \frac{b}{a}\]

2.4 Normalisation

Even in clean synthetic data set \(\rightarrow\) marginal precursor intensity distributions across samples not well aligned.

\(\rightarrow\) Can cause power loss & bias

2.4.1 Median centering

• Is a popular choice
• Often works well in very clean datasets
• But more advance normalisation often needed in real data.

Spike-in

Heart Example: Median Normalisation

• PXD006675 on PRIDE: Doll S, Dre, B, en M, Geyer PE, Itzhak DN, Braun C, Doppler SA, Meier F, Deutsch MA, Lahm H, Lange R, Krane M, Mann M, Region and cell-type resolved quantitative proteomic map of the human heart. Nat Commun, 8(1):1469(2017)

• 3 patients

• proteomes from biopsies of atrium and ventriculum in the left and the right heart region

• DDA

Why?

Median Normalisation

2.4.2 Median of Ratios method

• Borrowed from DESeq2: originally developed for RNA-seq data analysis
• Also correct for differences in composition of the proteomes in the different samples.

The method:

1. Calculates pseudo-reference sample (row-wise mean of log2 intensities, which corresponds to the log2 transformed geometric mean).
2. Calculates log2 ratios of each sample w.r.t. the pseudo-reference.
3. Column wise median of these log2 ratios is taken as the sample based normalisation factor on the log2 scale.

Figure 1: Illustration of the median-of-ratios normalisation on a toy example with four features and three samples. Starting from the matrix of log2 intensities (left), a pseudo-reference sample is obtained by taking the row-wise (feature-wise) mean. Log2 fold changes of every sample relative to this reference are then calculated (middle). Finally, the sample-specific normalisation factor is obtained as the median of these log2 fold changes over all features (right).

Spike-in

Heart Example: Median of Ratio Normalisation

compared to median normalisation:

2.5 Summarisation

• Summarise/Aggregate precursor-level intensities into a protein expression value.

• Peptide-level data for one of the spiked UPS proteins.

1. Protein can consist of many precursors
2. Huge effect of precursor species
3. Effect of spike-in concentration
4. Context dependent missingness
5. Often outliers, e.g. misidentification or fluctuations during MS acquisition.

2.5.1 Median Summarisation

• Problematic
• Biased FC \(\rightarrow\) Does not account for precursor species. Missingness.

2.5.2 MaxLFQ Summarisation

• Summarization should account for precursor species!
• Popular choice: maxLFQ

1. Calculates all pairwise log ratio’s between samples only using their shared precursors. Particularly, it uses the median of the log ratio’s between the shared precursors s, i.e.  \[r_{ij} = median(y_{sj}-y_{si})\].

2. Estimates the summaries by solving

\[ \sum_i\sum_j(y^{prot}_j - y^{prot}_i - r_{ij})^2 \]

Why?

3 Statistical Modeling and Inference

3.1 Motivating dataset

Movie clip

• Pathogen: causes tularemia
• Metabolic adaptation key for intracellular life cycle of pathogenic microorganisms.
• Upon entry into host cells quick phasomal escape and active multiplication in cytosolic compartment.
• Franciscella is auxotroph for several amino acids, including arginine.
• Inactivation of arginine transporter delayed bacterial phagosomal escape and intracellular multiplication.
• Experiment to assess difference in proteome using 3 WT vs 3 ArgP KO mutants

3.2 Import the data in R

Click to see code

1. Load libraries

library(tidyverse)
library(limma)
library(QFeatures)
library(msqrob2)
library(plotly)
library(ggplot2)

2. We use a peptides.txt file from MS-data quantified with maxquant that contains MS1 intensities summarized at the peptide level.

peptidesFile <- "https://raw.githubusercontent.com/statOmics/PDA22GTPB/data/quantification/francisella/peptides.txt"

3. Maxquant stores the intensity data for the different samples in columnns that start with Intensity. We can retreive the column names with the intensity data with the code below:

ecols <- grep("Intensity\\.", names(read.delim(peptidesFile)))

4. Read the data and store it in QFeatures object

qf <- readQFeatures(
  assayData = read.delim(peptidesFile),
  fnames = 1,
  quantCols =  ecols,
  name = "peptideRaw")

5. Update data with information on design

colData(qf)$genotype <- qf[[1]] |> 
  colnames() |> 
  substr(12,13) |>
  as.factor() |> 
  relevel("WT")
colData(qf)

DataFrame with 6 rows and 1 column
                         genotype
                         <factor>
Intensity.1WT_20_2h_n3_3       WT
Intensity.1WT_20_2h_n4_3       WT
Intensity.1WT_20_2h_n5_3       WT
Intensity.3D8_20_2h_n3_3       D8
Intensity.3D8_20_2h_n4_3       D8
Intensity.3D8_20_2h_n5_3       D8

3.3 Preprocessing

Click to see code to log-transfrom the data

1. Log transform

• peptides with zero intensities are missing peptides and should be represent with a NA value rather than 0.

qf <- zeroIsNA(qf, "peptideRaw") # convert 0 to NA

• Logtransform data with base 2

qf <- logTransform(qf, base = 2, i = "peptideRaw", name = "peptides_log")

2. Filtering

Only use proteotypic proteins

qf <- filterFeatures(qf, ~ !grepl(pattern = ";", Proteins))

• Remove reverse sequences (decoys) and contaminants. Note that this is indicated by the column names Reverse and depending on the version of maxQuant with Potential.contaminants or Contaminants.

qf <- filterFeatures(qf,~Reverse != "+")
qf <- filterFeatures(qf,~ Contaminant != "+")

• We keep peptides that were observed at last 4 times out of the \(n = 6\) samples, so that we can estimate the peptide characteristics.

We tolerate the following proportion of NAs: \(\text{pNA} = \frac{(n - 4)}{n} = 0.33\), so we keep peptides that are observed in at least 66% of the samples. This is an arbitrary value that may need to be adjusted depending on the experiment and the data set.

nObs <- 4
n <- ncol(qf[["peptides_log"]])

(qf <- filterNA(qf, i = "peptides_log", pNA = (n - nObs) / n))

An instance of class QFeatures (type: bulk) with 2 sets:

 [1] peptideRaw: SummarizedExperiment with 10461 rows and 6 columns 
 [2] peptides_log: SummarizedExperiment with 5052 rows and 6 columns 

We keep 5052 peptides upon filtering.

3. Normalization by median centering

qf <-  sweep( #4. Subtract log2 norm factor column-by-column (MARGIN = 2)
  qf,
  MARGIN = 2, 
  STATS = nfLogMedianOfRatios(qf, i="peptides_log"), #1.
  i = "peptides_log",
  name = "peptides_norm"
)

4. Summarization. We use maxLFQ from the iq package.

qf <- aggregateFeatures(qf,
    i = "peptides_norm", 
    fcol = "Proteins", 
    name = "proteins",
    fun = function(X) iq::maxLFQ(X)$estimate
)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

Plot of preprocessed data

qf[, , "peptides_norm"] |> #1.
  longForm(colvars = c( "genotype")) |>  #2. 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + #3.
  aes(x = value,
      colour = genotype,
      group = colname) +
  geom_density() +
  theme_minimal() +
  ggtitle("Peptide-level")

qf[, , "proteins"] |> #1.
  longForm(colvars = c( "genotype")) |>  #2. 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + #3.
  aes(x = value,
      colour = genotype,
      group = colname) +
  geom_density() +
  theme_minimal() +
    ggtitle("Protein-level")

3.4 Summarized data structure

3.4.1 Design

qf |> 
  colData() |> 
  knitr::kable()

	genotype
Intensity.1WT_20_2h_n3_3	WT
Intensity.1WT_20_2h_n4_3	WT
Intensity.1WT_20_2h_n5_3	WT
Intensity.3D8_20_2h_n3_3	D8
Intensity.3D8_20_2h_n4_3	D8
Intensity.3D8_20_2h_n5_3	D8

• WT vs KO
• 3 vs 3 repeats

3.4.2 Summarized intensity matrix

qf[["proteins"]] |> 
  assay() |> 
  head() |> 
  knitr::kable()

	Intensity.1WT_20_2h_n3_3	Intensity.1WT_20_2h_n4_3	Intensity.1WT_20_2h_n5_3	Intensity.3D8_20_2h_n3_3	Intensity.3D8_20_2h_n4_3	Intensity.3D8_20_2h_n5_3
WP_003013731	25.91971	26.03992	26.27133	25.93598	26.19855	26.22295
WP_003013909	25.87033	25.65042	25.74652	26.07970	25.90919	26.03624
WP_003014068	26.79067	26.76619	27.08102	26.62292	26.80336	26.77253
WP_003014122	25.30922	24.96019	25.18033	24.92020	24.80622	24.86137
WP_003014123	25.86714	25.66709	25.77431	25.70902	25.71849	25.53778
WP_003014302	27.85088	27.93109	27.87884	28.01439	28.02191	28.05122

• 1020 proteins

3.5 Hypothesis testing: a single protein

Movie clip

3.5.1 T-test

\[ \text{Effect size:} \log_2 \text{FC} = \bar{y}_{p1}-\bar{y}_{p2} \]

\[ T_g=\frac{\log_2 \widehat{\text{FC}}}{\text{se}_{\log_2 \widehat{\text{FC}}}} \]

\[ T_g=\frac{\widehat{\text{signal}}}{\widehat{\text{Noise}}} \]

If we can assume equal variance in both treatment groups:

\[ \text{se}_{\log_2 \text{FC}}=\text{SD}\sqrt{\frac{1}{n_1}+\frac{1}{n_2}} \]

3.5.2 Linear model

\[ \begin{array}{lcl} y_i &=& \beta_0 + \beta_{D8}x_{i,D8} + \epsilon_i\\ \epsilon_i &\sim& N(0,\sigma^2)\\ \end{array} \] - \(x_{i,D8}\) is dummy variable

\[ \begin{array}{lcl} x_{i,D8} &=& 0 \text{ (WT)}\\ x_{i,D8} &=& 1 \text{ (D8)} \end{array} \] - \(\log_2\) FC?

\[ \begin{array}{lcl} \hat \mu_{WT} &=& \hat\beta_0 \\ \hat \mu_{D8} &=& \hat\beta_0 + \hat{\beta}_{D8}\\ \log_2 \widehat{\text{FC}} &=& \hat{\beta}_{D8} \end{array} \]

In R

lm1  <- lm(intensity ~ genotype, data = WP_003023392)
summary(lm1)

Call:
lm(formula = intensity ~ genotype, data = WP_003023392)

Residuals:
       1        2        3        4        5        6 
 0.07266 -0.12677  0.05411  0.03622 -0.15565  0.11943 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 26.10746    0.07308  357.27 3.68e-10 ***
genotypeD8  -1.36902    0.10334  -13.25 0.000188 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1266 on 4 degrees of freedom
Multiple R-squared:  0.9777,    Adjusted R-squared:  0.9721 
F-statistic: 175.5 on 1 and 4 DF,  p-value: 0.0001877

\[ t=\frac{\log_2\widehat{\text{FC}}}{\text{se}_{\log_2\widehat{\text{FC}}}}=\frac{-1.37}{0.103}=-13.2 \]

• Is t = -13.2 indicating that there is an effect?

• How likely is it to observe t = -13.2 when there is no effect of the argP KO on the protein expression?

3.5.3 Null hypothesis (\(H_0\)) and alternative hypothesis (\(H_1\))

• With data we can never prove a hypothesis (falsification principle of Popper)

• With data we can only reject a hypothesis

• In general we start from alternative hypothese \(H_1\): we want to show an effect of the KO on abundance of a protein

\(H_1\): On average the protein abundance in WT is different from that in KO

• But, we will assess this by falsifying the opposite: \(H_0\): On average the protein abundance in WT is equal to that in KO<-

lm1 <- lm(intensity ~ genotype, data = WP_003023392)
summary(lm1)

Call:
lm(formula = intensity ~ genotype, data = WP_003023392)

Residuals:
       1        2        3        4        5        6 
 0.07266 -0.12677  0.05411  0.03622 -0.15565  0.11943 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 26.10746    0.07308  357.27 3.68e-10 ***
genotypeD8  -1.36902    0.10334  -13.25 0.000188 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1266 on 4 degrees of freedom
Multiple R-squared:  0.9777,    Adjusted R-squared:  0.9721 
F-statistic: 175.5 on 1 and 4 DF,  p-value: 0.0001877

• How likely is it to observe an equal or more extreme effect than the one observed in the sample when the null hypothesis is true?

• When we make assumptions about the distribution of our test statistic we can quantify this probability: p-value. The p-value will only be calculated correctly if the underlying assumptions hold!

• When we repeat the experiment, the probability to observe a fold change for this protein that is more extreme than a 2.58 fold (\(\log_2 FC=-1.37\)) down or up regulation by random change (if \(H_0\) is true) is 2 out of 10000.
  

• If the p-value is below a significance threshold \(\alpha\) we reject the null hypothesis. We control the probability on a false positive result at the \(\alpha\)-level (type I error)

• Note, that the p-values are uniform under the null hypothesis, i.e. when \(H_0\) is true all p-values are equally likely.

3.6 Multiple hypothesis testing

Movie clip

• Consider testing DA for all \(m=1020\) proteins simultaneously

• What if we assess each individual test at level \(\alpha\)? \(\rightarrow\) Probability to have a false positive (FP) among all m simultatenous test \(>>> \alpha= 0.05\)

• Indeed for each non DA protein we have a probability of 5% to return a FP.

• In a typical experiment the majority of the proteins are non DA.

• So an upperbound of the expected FP is \(m \times \alpha\) or \(1020 \times 0.05=51\).

\(\rightarrow\) Hence, we are bound to call many false positive proteins each time we run the experiment.

3.6.1 Multiple testing: False discovery rate

Movie clip

• FDR: Expected proportion of false positives on the total number of positives you return.
• An FDR of 1% means that on average we expect 1% false positive proteins in the list of proteins that are called significant.
• Defined by Benjamini and Hochberg in their seminal paper Benjamini, Y. and Hochberg, Y. (1995). “Controlling the false discovery rate: a practical and powerful approach to multiple testing”. Journal of the Royal Statistical Society Series B, 57 (1): 289–300.

The False Discovery Proportion (FDP) is the fraction of false positives that are returned, i.e. 

\[ FDP = \frac{FP}{R} \]

• However, this quantity cannot be observed because in practice we only know the number of proteins for which we rejected \(H_0\), \(R\).
• But, we do not know the number of false positives, \(FP\).

Therefore, Benjamini and Hochberg, 1995, defined The False Discovery Rate (FDR) as \[ \text{FDR} = \text{E}\left[\frac{FP}{R}\right] =\text{E}\left[\text{FDP}\right] \] the expected FDP.

• Controlling the FDR allows for more discoveries (i.e. longer lists with significant results), while the fraction of false discoveries among the significant results in well controlled on average. As a consequence, more of the true positive hypotheses will be detected.

Intuition of BH-FDR procedure

Consider \(m = 1000\) tests

• Suppose that a researcher rejects all null hypotheses for which \(p < 0.01\).

• If we use \(p < 0.01\), we expect \(0.01 \times m_0\) tests to return false positives.

• A conservative estimate of the number of false positives that we can expect can be obtained by considering that the null hypotheses are true for all features, \(m_0 = m = 1000\).

• We then would expect \(0.01 \times 1000 = 10\) false positives (\(FP=10\)).

• Suppose that the researcher found 200 proteins with \(p<0.01\) (\(R=200\)).

• The proportion of false positive results (FDP = false positive proportion) among the list of \(R=200\) proteins can then be estimated as

\[ \widehat{\text{FDP}}=\frac{\widehat{FP}}{R}=\frac{10}{200}=\frac{0.01 \times 1000}{200} = 0.05. \]

Benjamini and Hochberg (1995) procedure for controlling the FDR at \(\alpha\)

1. Let \(p_{(1)}\leq \ldots \leq p_{(m)}\) denote the ordered \(p\)-values.

2. Find the largest integer \(k\) so that \[ \frac{p_{(k)} \times m}{k} \leq \alpha \] \[\text{or}\] \[ p_{(k)} \leq k \times \alpha/m \]

3. If such a \(k\) exists, reject the \(k\) null hypotheses associated with \(p_{(1)}, \ldots, p_{(k)}\). Otherwise none of the null hypotheses is rejected.

The adjusted \(p\)-value (also known as the \(q\)-value in FDR literature): \[ q_{(i)}=\tilde{p}_{(i)} = \min\left[\min_{j=i,\ldots, m}\left(m p_{(j)}/j\right), 1 \right]. \] In the hypothetical example above: \(k=200\), \(p_{(k)}=0.01\), \(m=1000\) and \(\alpha=0.05\).

Francisella Example

Movie clip

Click to see code

ttestMx <- function(y,group) {
    test <- try(t.test(y[group],y[!group],var.equal=TRUE),silent=TRUE)
    if(is(test,"try-error")) {
      return(c(log2FC=NA,se=NA,tstat=NA,p=NA))
      } else {
      return(c(log2FC= (test$estimate%*%c(1,-1)),se=test$stderr,tstat=test$statistic,pval=test$p.value))
      }
 }
 
 res <- apply(
    assay(qf[["proteins"]]), 
    1, 
    ttestMx,
    group = colData(qf)$genotype=="D8") |> 
  t()
 colnames(res) <- c("logFC","se","tstat","pval")
 res <- res |> 
   as.data.frame() |> 
   na.exclude() |> 
   arrange(pval)
 res$adjPval <- p.adjust(res$pval, "fdr")
 alpha <- 0.05
res$adjAlphaForm <- paste0(1:nrow(res)," x ",alpha,"/",nrow(res))
res$adjAlpha <- alpha * (1:nrow(res))/nrow(res) 
res$"pval < adjAlpha" <- res$pval < res$adjAlpha 
res$"adjPval < alpha" <- res$adjPval < alpha 

	logFC	pval	adjPval	adjAlphaForm	adjAlpha	pval < adjAlpha	adjPval < alpha
WP_003034781	0.4224142	0.0000022	0.0022721	1 x 0.05/1020	0.0000490	TRUE	TRUE
WP_003040849	-1.2263612	0.0000076	0.0026390	2 x 0.05/1020	0.0000980	TRUE	TRUE
WP_011733723	-0.3258931	0.0000078	0.0026390	3 x 0.05/1020	0.0001471	TRUE	TRUE
WP_003041244	0.2057070	0.0000791	0.0201736	4 x 0.05/1020	0.0001961	TRUE	TRUE
WP_003038430	-0.3836478	0.0001838	0.0280097	5 x 0.05/1020	0.0002451	TRUE	TRUE
WP_003023392	-1.3690155	0.0001877	0.0280097	6 x 0.05/1020	0.0002941	TRUE	TRUE
WP_011733588	-0.3823744	0.0001922	0.0280097	7 x 0.05/1020	0.0003431	TRUE	TRUE
WP_011733645	0.3040827	0.0002362	0.0301195	8 x 0.05/1020	0.0003922	TRUE	TRUE
WP_004339068	0.1975424	0.0006535	0.0730468	9 x 0.05/1020	0.0004412	FALSE	FALSE
WP_003018840	-0.5013825	0.0007161	0.0730468	10 x 0.05/1020	0.0004902	FALSE	FALSE

...

WP_003022463	-0.0003682	0.9989893	0.9999324	1019 x 0.05/1020	0.049951	FALSE	FALSE
WP_003033335	-0.0000299	0.9999324	0.9999324	1020 x 0.05/1020	0.050000	FALSE	FALSE

Results

Click to see code

volcanoT <- plotVolcano(res) 

volcanoT

3.7 Moderated Statistics

Movie clip

Problems with ordinary t-test

Click to see code

problemPlots <- list() 
problemPlots[[1]] <- res |> 
  ggplot(aes(x = logFC, y = se, color = adjPval < 0.05)) +
    geom_point(cex = 2.5) +
    scale_color_manual(values = alpha(c("black", "red"), 0.5)) +
    theme_minimal() 

for (i in 2:3)
{
problemPlots[[i]] <- colData(qf) |> 
    as.data.frame() |> 
    mutate(intensity = qf[["proteins"]][rownames(res)[i],] |> 
             assay() |> 
             c()) |> 
    ggplot(aes(x=genotype,y=intensity)) +
    geom_point() + 
    ggtitle(rownames(res)[i]) +
    ylim(c(25.8,27.4)) +
      theme_minimal() 
}

problemPlots

[[1]]

[[2]]

[[3]]

3.7.1 Empirical Bayes

Figure courtesy to Rafael Irizarry

\[ T_p^{mod} = \frac{\log_2 \widehat{\text{FC}} }{C \quad \tilde{S}_p} , \]

• empirical Bayes theory provides formal framework for borrowing strength across proteins,
• Implemented in popular bioconductor package limma and msqrob2

\[ \tilde{S}_p=\sqrt{\frac{d_pS_p^2+d_0S_0^2}{d_p+d_0}}, \]

• \(S_0^2\): common variance (over all proteins)
• Moderated t-statistic is t-distributed with \(d_0+d_p\) degrees of freedom.
• Note that the degrees of freedom increase by borrowing strength across proteins!

Click to see the code

1. We model the protein level expression values using the msqrob function. By default msqrob2 estimates the model parameters using robust regression.

We will model the data with a different group mean for every genotype. The group is incoded in the variable genotype of the colData. We can specify this model by using a formula with the factor genotype as its predictor: formula = ~genotype.

Note, that a formula always starts with a symbol ‘~’.

qf <- msqrob(object = qf, i = "proteins", formula = ~genotype)

2. Inference

We first explore the design of the model that we specified using the the package ExploreModelMatrix

library(ExploreModelMatrix)
VisualizeDesign(colData(qf),~genotype)$plotlist[[1]]

We have two model parameters, the (Intercept) and genotypeD8. This results in a model with two group means:

1. For the wild type (WT) the expected value (mean) of the log2 transformed intensity y for a protein will be modelled using

\[\text{E}[Y\vert \text{genotype}=\text{WT}] = \text{(Intercept)}\]

2. For the knockout genotype D8 the expected value (mean) of the log2 transformed intensity y for a protein will be modelled using

\[\text{E}[Y\vert \text{genotype}=\text{D8}] = \text{(Intercept)} + \text{genotypeD8}\]

The average log2FC between D8 and WT is thus \[\log_2\text{FC}_{D8-WT}= \text{E}[Y\vert \text{genotype}=\text{D8}] - \text{E}[Y\vert \text{genotype}=\text{WT}] = \text{genotypeD8} \]

Hence, assessing the null hypothesis that there is no differential abundance between D8 and WT can be reformulated as

\[H_0: \text{genotypeD8}=0\] We can implement a hypothesis test for each protein in msqrob2 using the code below:

L <- makeContrast("genotypeD8 = 0", parameterNames = c("genotypeD8"))
qf <- hypothesisTest(object = qf, i = "proteins", contrast = L)

We can show the list with all significant DE proteins at the 5% FDR using

msqrobCollect(qf[["proteins"]],L) |> 
  arrange(pval) |>
  filter(adjPval<0.05)

                             logFC         se       df          t         pval
genotypeD8.WP_003040849 -1.2263612 0.08020151 6.087933 -15.290999 4.358069e-06
genotypeD8.WP_003023392 -1.3690155 0.11105902 6.087933 -12.326919 1.560512e-05
genotypeD8.WP_003033719 -1.1615474 0.12602389 5.828253  -9.216882 1.087304e-04
genotypeD8.WP_003026016 -0.9540456 0.11530658 6.080135  -8.273991 1.573910e-04
genotypeD8.WP_003033046 -0.8219170 0.10214460 5.815009  -8.046603 2.307079e-04
genotypeD8.WP_003038350  0.8077325 0.10952806 6.087933   7.374663 2.978545e-04
                            adjPval   contrast      feature
genotypeD8.WP_003040849 0.004301414 genotypeD8 WP_003040849
genotypeD8.WP_003023392 0.007701128 genotypeD8 WP_003023392
genotypeD8.WP_003033719 0.035772317 genotypeD8 WP_003033719
genotypeD8.WP_003026016 0.038836240 genotypeD8 WP_003026016
genotypeD8.WP_003033046 0.045541744 genotypeD8 WP_003033046
genotypeD8.WP_003038350 0.048997070 genotypeD8 WP_003038350

We can also visualise the results using a volcanoplot

volcano <- plotVolcano(
    rowData(qf[["proteins"]])$genotypeD8) +
  ggtitle("msqrob2")

• The volcano plot opens up when using the EB variance estimator

• Borrowing strength to estimate the variance using empirical Bayes solves the issue of returning proteins with a low fold change as significant due to a low variance.

3.7.2 Shrinkage of the variance and moderated t-statistics

• Small variances are shrunken towards the common variance resulting in larger EB variance estimates
• Large variances are shrunken towards the common variance resulting in smaller EB variance estimates
• Pooled degrees of freedom of the EB variance estimator are larger because information is borrowed across proteins to estimate the variance

3.8 Plots

4 Some Concepts on Experimental Design

4.1 Random Sampling

• Random sampling is closely related to the concept of the population or the scope of the study.

• Based on a sample of subjects, the researchers want to come to conclusions that hold for

  • all kinds of people
  • people with hypertension

• Scope of the study should be well specified before the start of the study.

• Representative sample: For the statistical analysis to be valid, it is required that the subjects are selected completely at random from the population to which we want to generalize our conclusions.

• Selecting completely at random from a population implies:

  • all subjects in the population should have the same probability of being selected in the sample,
  • the selection of a subject in the sample should be independent from the selection of the other subjects in the sample.

4.2 Randomisation

• Make sure that groups are comparable / Avoid systematic differences between groups
• Randomisation: treatments of interest are attributed at random to the experimental units (subjects, animals, plants, cages, …)

4.3 Control

Good control is nessecary

\(\rightarrow\) Placebo controlled double blind experiments
\(\rightarrow\) Injection of control animal with blank containing same solvents, etc.

4.4 Replication / Sample size

Within an experiment: Enables to estimate uncertainty / biological variability

• Replication is essential for quantifying the noise

• Noise: biological and technical in nature

Between experiments: Any true finding should be reproducible

• Genuine replicates include all sources of variability: technical + biological

• Technical replicates are important if you assess new technologies

Movie clip

\[ \log_2 \text{FC} = \bar{y}_{p1}-\bar{y}_{p2} \]

\[ T_g=\frac{\log_2 \text{FC}}{\text{se}_{\log_2 \text{FC}}} \]

\[ T_g=\frac{\widehat{\text{signal}}}{\widehat{\text{Noise}}} \]

If we can assume equal variance in both treatment groups:

\[ \text{se}_{\log_2 \text{FC}}=\text{SD}\sqrt{\frac{1}{n_1}+\frac{1}{n_2}} \]

\(\rightarrow\) if number of bio-repeats increases we have a higher power!

• cfr. Study of tamoxifen treated Estrogen Recepter (ER) positive breast cancer patients (tutorial)

4.5 Randomized complete block designs

\[\sigma^2= \sigma^2_{bio}+\sigma^2_\text{lab} +\sigma^2_\text{extraction} + \sigma^2_\text{run} + \ldots\]

• Biological: fluctuations in protein level between mice, fluctations in protein level between cells, …
• Technical: cage effect, lab effect, week effect, plasma extraction, MS-run, …

4.6 Nature methods: Points of significance - Blocking

https://www.nature.com/articles/nmeth.3005.pdf

4.6.1 Mouse example

Duguet et al. (2017) MCP 16(8):1416-1432. doi: 10.1074/mcp.m116.062745

Subsampling to show impact of design

Completely randomised design (CRD):

• 4 mice vs 4 mice
• Only one cell type per mouse
• Linear model:

\[ \sim \text{type} \]

Randomised complete block design:

• 4 mice
• Both regulatory and ordinary T-cells per mouse

• All treatments of interest are present within block!
• We can estimate the effect of the treatment within block!
• We can isolate the between block variability from the analysis using linear model: \[ \sim \text{type} + \text{mouse} \]

RCB design (blocking):

Results

• Advantage if block effect/variability between mice is large. Big reduction in residual errors.
  
• Otherwise loss in degrees of freedom (more model parameters have to be estimated).

\[ \text{SD}_\text{res} = \sqrt{\frac{\sum e_i^2}{n-p}} \]

with \(n-p\) equals

• 8-5 = 3 for RCB and
• 8 - 2 = 6 for CRD.

5 Take home messages

• Explore your data in every stage of the analysis

• Data processing:

  • Avoid imputation
  • Filtering
  • log2 transformation: mean-variance + log FC
  • More advanced normalisation methods that can account for compositional differences are often needed for real data
  • Summarisation should account for differences in precursor/peptide characteristics

• Statistical modeling and inference:

  • Carefully think about sources of variation in experiment
  • Translate research question in contrast/effect size
  • Multiple testing
  • Empirical Bayes: stabilise variance estimation
  • Critically explore results

• Experimental design:

  • Thinking about experimental design and data analysis prior to the experiment is key!

  • Random sampling

  • Randomisation

  • Good control

  • Sample size

  • Blocking: very powerful approach to isolate important sources of variability from the experiment.

    • Both biological, e.g. animal, cell culture
    • Technical, e.g. batch effects, …

6 References

Staes, An, Teresa Mendes Maia, Sara Dufour, et al. 2024. “Benefit of in Silico Predicted Spectral Libraries in Data‑independent Acquisition Data Analysis Workflows.” Journal of Proteome Research 23 (6): 2078–89. https://doi.org/10.1021/acs.jproteome.4c00048.