limma_proteomics_tutorial_2025.Rmd

---
title: "Limma for Proteomics"
author: "Blythe Durbin-Johnson, Ph.D."
output:
  html_document:
    keep_md: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Introduction

Limma is an R package (developed for use with gene expression microarrays) that is used for differential abundance/expression analysis of proteomics, metabolomics, RNA sequencing, and other 'omics data.

### Core steps of limma analysis

1. A linear model (e.g. ANOVA or regression) is fitted to each protein.

2. Empirical Bayes smoothing is used to get better estimates of standard errors of log fold changes, which are then used in differential abundance testing.

3. P-values are adjusted for multiple testing across proteins using the Benjamini-Hochberg false discovery rate controlling method.

Limma allows for more complicated experimental designs and models than many other packages.

References:

Matthew E. Ritchie, Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, Gordon K. Smyth, limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Research, Volume 43, Issue 7, 20 April 2015, Page e47, https://doi.org/10.1093/nar/gkv007

Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3. doi: 10.2202/1544-6115.1027. Epub 2004 Feb 12. PMID: 16646809.

https://www.bioconductor.org/packages//2.11/bioc/vignettes/limma/inst/doc/usersguide.pdf

### Data

We will use proteomics data from Kistner et al. on 50 subjects with chronic lymphocytic leukemia (CLL).  

The differential abundance analysis will compare subjects with mutated (M) IGHV versus unmutated (U) IGHV.  CLL patients with mutated IGHV tend to have a better prognosis.


References:

Franziska Kistner, Justus L. Grossmann, Ludwig R. Sinn, Vadim Demichev
bioRxiv 2023.06.20.545604; doi: https://doi.org/10.1101/2023.06.20.545604

## Analysis

Load R packages:

```{r, warning = FALSE, message = FALSE}
library(dplyr)
library(limma)
library(ggplot2)
library(ggrepel)
library(ComplexHeatmap)
```

### Read in proteomics data and metadata

Read in proteomics data, clean up:
```{r}
dat0 <- read.delim("legacy-report.pg_matrix.tsv", check.names = FALSE)

# Split out protein annotation and protein abundance data
anno <- select(dat0, Protein.Group:First.Protein.Description)
dat <- select(dat0, -(Protein.Group:First.Protein.Description))

# Clean up column names so they match filenames in metadata
colnames(dat) <- gsub("D:\\PXD022216\\", "", colnames(dat), fixed = TRUE)

# Set rownames to Protein.Group
rownames(dat) <- anno$Protein.Group
```

Read in metadata, reorder so matches order of proteomics data.  *Very important!*

```{r}
pdata <- read.delim("Annotation_file_timsTOF.txt")
rownames(pdata) <- pdata$File.name

# Reorder to match proteomics data
pdata <- pdata[colnames(dat),]

# Check again to make sure metadata matches proteomics data
identical(pdata$File.name, colnames(dat))
```
Set missing values to 0:

```{r}
na.mat <- is.na(dat) 
dat[na.mat] <- 0
```

We will log transform the proteomics data prior to analysis.  This helps make the data closer to normally distributed and makes the variability more constant between low and high expressed proteins.

Log transform data:

```{r}
dat.trans <- log2(dat + 1)
```

Boxplot of data to evaluate need for between-sample normalization:

```{r}
boxplot(dat.trans, xaxt = "n")
```

Between sample normalization appears to have already been performed.  We could use the limma function normalizeCyclicLoess to normalize if this wasn't the case:

Example of normalization (not run):

```{r, eval = FALSE}
dat.trans <- normalizeCyclicLoess(dat.trans)
```

Next, we filter the data to remove proteins present in few samples.  This helps with the following:

1. Testing fewer proteins makes the multiple testing adjustment less strict (the remaining proteins have a better chance of being significant).

2. Having proteins that are absent (zero) in lots of samples introduces large numbers of ties in the data, causing the variability to be underestimated.

Filter to proteins present in at least half of samples:

```{r}
num.nonmissing <- rowSums(!na.mat)
keep <- num.nonmissing >= ncol(dat.trans)/2
dat.filtered <- dat.trans[keep,]
```

### Multidimensional scaling plots

A multidimensional scaling (MDS) plot is a two dimensional view of the data showing the relative distances between protein profiles.  

We will use the limma function plotMDS to get the MDS coordinates then make the actual plot in ggplot2 (https://ggplot2.tidyverse.org/) for extra flexibility.


```{r}
coords <- plotMDS(dat.filtered, plot = FALSE)
plotdat <- pdata
plotdat$MDS_1 <- coords$x
plotdat$MDS_2 <- coords$y
```

#### MDS plot by IGHV status
```{r}
ggplot(plotdat, aes(x = MDS_1, y = MDS_2, color = IGHV.status)) + geom_point(size = 3) + theme_bw() + labs(color = "IGHV status")
```

#### MDS plot by gender
```{r}
ggplot(plotdat, aes(x = MDS_1, y = MDS_2, color = Gender)) + geom_point(size = 3) + theme_bw()
```

#### MDS plot by leukocyte count
```{r}
ggplot(plotdat, aes(x = MDS_1, y = MDS_2, color = Leukocyte.count)) + geom_point(size = 3) + theme_bw() + scale_color_viridis_c() + labs(color = "Leukocyte count")
```

### Differential abundance analysis in limma

We will compare expression of each protein between subjects with mutated and unmutated IGHV.

The design matrix is how we specify the experimental design in limma.

Set up design matrix:

```{r}
mm <- model.matrix(~0 + IGHV.status, data = pdata)
```

Fit model to each protein:

```{r}
fit <- lmFit(dat.filtered, mm)
```

Specify comparison ("contrast") being tested:

```{r}
contr <- makeContrasts(IGHV.statusM - IGHV.statusU, levels = colnames(coef(fit)))
contr2 <- contrasts.fit(fit, contrasts = contr)
```

Empirical Bayes smoothing adjusts the standard error of each protein to be closer to the average standard error, compensating for unusually high or low variability proteins.

Empirical Bayes smoothing:

```{r}
contr3 <- eBayes(contr2)
```

Finally, we adjust for testing multiple proteins.  With `r nrow(dat.filtered)` proteins after filtering, we would expect `r 0.05*nrow(dat.filtered)` to have a raw p-value less than 0.05, even if there were no true differences.

Defining significant as FDR adjusted P < 0.05 means that we expect 5% of the proteins with an adjusted p-value < 0.05 to be false discoveries.

Display differential abundance results with Benjamini-Hochberg false discovery rate adjusted p-values:

```{r}
results <- topTable(contr3, sort.by = "P", n = Inf)
head(results)
```

Results include the following columns:

* logFC: log2 fold change for IGHV status M - U

* AveExpr: Average protein level across all samples, on transformed scale

* t: t statistic, or log fold change divided by its standard error

* P.Value: Raw p-value from the test that the logFC differs from 0

* adj.P.Val: Benjamini-Hochberg false discovery rate adjusted p-value

* B: log odds of differential abundance


How many proteins differ significantly between mutated and unmutated IGHV status subjects?:

```{r}
length(which(results$adj.P.Val < 0.05))
```

Merge in annotation information:

```{r}
results$Protein.Group <- rownames(results)
results <- left_join(results, anno, by = "Protein.Group")
head(results)
```

### Volcano plot of results
A volcano plot shows -log10(p-value) plotted against the log fold change.

```{r}
results$isSig <- ifelse(results$adj.P.Val < 0.05, "Significant", "Not Significant")
top.results <- results[1:15,] # we will only label the most significant proteins
ggplot(results, aes(x = logFC, y = -log10(P.Value), col = isSig)) + geom_point() + 
  scale_color_manual(values = c("grey", "blue")) + 
  geom_text_repel(data = top.results, aes(label = Genes), max.overlaps = Inf, color = "black") + 
  theme_bw() + 
  labs(color = NULL)
```


### Heatmap of top 20 proteins

A heatmap displays the sample-level expression of a set of proteins.

We'll use the R package ComplexHeatmap (https://jokergoo.github.io/ComplexHeatmap-reference/book/)


```{r}
toplot <- results$Protein.Group[1:20] 
plotdat <- as.matrix(dat.filtered[toplot,])
new.rownames <- anno$Genes[match(rownames(plotdat), anno$Protein.Group)]

set.seed(99) # Annotation bar colors are random!
col_ha <- HeatmapAnnotation(`IGHV status` = pdata$IGHV.status, annotation_name_side = "left")

Heatmap(plotdat, heatmap_legend_param = list(title = "log2\nNormalized\nAbundance"), show_column_names = FALSE, row_labels = new.rownames, bottom_annotation = col_ha)
```

Plot using Z scores:

```{r}
toplot <- results$Protein.Group[1:20] 
plotdat <- as.matrix(dat.filtered[toplot,])
# convert to Z scores
plotdat <- t(apply(plotdat, 1, scale))

new.rownames <- anno$Genes[match(rownames(plotdat), anno$Protein.Group)]

set.seed(99) # Annotation bar colors are random!
col_ha <- HeatmapAnnotation(`IGHV status` = pdata$IGHV.status, annotation_name_side = "left")

Heatmap(plotdat, heatmap_legend_param = list(title = "Z score"), show_column_names = FALSE, row_labels = new.rownames, bottom_annotation = col_ha)
```

### Boxplots of individual proteins
```{r}
plotdat <- pdata
plotdat$L1TD1 <- as.numeric(dat.filtered["Q5T7N2",])
plotdat$LMNA <- as.numeric(dat.filtered["P02545",])
```

```{r}
ggplot(plotdat, aes(x = IGHV.status, y = L1TD1, fill = IGHV.status)) + geom_boxplot() + theme_bw() + theme(legend.position = "none") + labs(x = "IGHV status")
```

```{r}
ggplot(plotdat, aes(x = IGHV.status, y = LMNA, fill = IGHV.status)) + geom_boxplot() + theme_bw() + theme(legend.position = "none") + labs(x = "IGHV status")
```

### R session information
```{r}
sessionInfo()
```