Background

Proteomics data were generated from the same samples as the transcriptomics data (3 time points × 6 replicates = 18 samples). The investigators provided log₂ transformed and filtered data. They also provided us with “raw” data, which has not yet been transformed but has been filtered. We will touch on the structure of the raw data and how you may want to prefilter and transform the data.

Similar to our transcriptomics analysis from Lesson 1, we will:

examine the structure of the data
add metadata
transform the data
transform data to fit assumptions for analysis
fit differential abundance model
apply cutoffs to identify differentially abundant proteins
visualize the results

Setup

Load packages

Load packages and define vector of custom colors.

library(dplyr)
library(limma)
library(ggplot2)
library(ggrepel)
library(ComplexHeatmap)

custom_colors <- c("salmon", "darkolivegreen3", "turquoise3", "mediumpurple1")

Define volcano_plot() function which we will use later to visualize differential expression results.

volcano_plot <- function(df, title) {

  df$negLogP <- -log10(df$adj.P.Val)

  fc_thresh <- 0.75
  p_thresh  <- -log10(0.01)

  df$Significant <- with(df, abs(logFC) >= fc_thresh & negLogP >= p_thresh)

  ggplot(df, aes(x = logFC, y = negLogP)) +
    geom_point(aes(color = Significant), alpha = 0.6) +
    
    geom_vline(xintercept = c(-fc_thresh, fc_thresh), 
               linetype = "dashed", color = "red") +
    geom_hline(yintercept = p_thresh, 
               linetype = "dashed", color = "blue") +
  
    geom_text_repel(
      data = subset(df, Significant),
      aes(label = Gene),
      size = 3,
      max.overlaps = 20,
      box.padding = 0.4,
      point.padding = 0.3
    ) +
    
    theme_minimal() +
    labs(title = title,
         x = "log2 Fold Change",
         y = "-log10(Adjusted p-value)") +
    scale_color_manual(values = c("FALSE" = "grey50", "TRUE" = "red"))

}

Import data

Create a character variable with the path to the tab-separated text file that contains the untransformed proteome data.

raw <- "data/Imputed_proteinGroups_perseus.txt"

If you’re running this code on your own machine, you will need to read the data from the GitHub Pages URL where it is hosted.

raw <- url("https://usda-ree-ars.github.io/SEAStatsData/omics2026/Imputed_proteinGroups_perseus.txt")

Read the untransformed data from the tab-separated text file.

data <- read.table(raw, 
                   sep = '\t', 
                   header = TRUE, 
                   comment.char = '#', 
                   quote = "", 
                   fill = TRUE, 
                   check.names = FALSE)

Filter and transform data

Filter the data

Let’s look at the column names. The first 18 columns of this dataframe are the intensities for each sample. There are a lot of columns that come from the preprocessing of the proteomics data.

colnames(data)

##  [1] "0h_01"                               
##  [2] "0h_02"                               
##  [3] "0h_03"                               
##  [4] "0h_04"                               
##  [5] "0h_05"                               
##  [6] "0h_06"                               
##  [7] "4h_01"                               
##  [8] "4h_02"                               
##  [9] "4h_03"                               
## [10] "4h_04"                               
## [11] "4h_05"                               
## [12] "4h_06"                               
## [13] "24h_01"                              
## [14] "24h_02"                              
## [15] "24h_03"                              
## [16] "24h_04"                              
## [17] "24h_05"                              
## [18] "24h_06"                              
## [19] "Only identified by site"             
## [20] "Reverse"                             
## [21] "Potential contaminant"               
## [22] "Peptides"                            
## [23] "Razor + unique peptides"             
## [24] "Unique peptides"                     
## [25] "Sequence coverage [%]"               
## [26] "Unique + razor sequence coverage [%]"
## [27] "Unique sequence coverage [%]"        
## [28] "Mol. weight [kDa]"                   
## [29] "Q-value"                             
## [30] "Score"                               
## [31] "Intensity"                           
## [32] "iBAQ peptides"                       
## [33] "MS/MS count"                         
## [34] "Peptides 0h_01"                      
## [35] "Peptides 0h_02"                      
## [36] "Peptides 0h_03"                      
## [37] "Peptides 0h_04"                      
## [38] "Peptides 0h_05"                      
## [39] "Peptides 0h_06"                      
## [40] "Peptides 4h_01"                      
## [41] "Peptides 4h_02"                      
## [42] "Peptides 4h_03"                      
## [43] "Peptides 4h_04"                      
## [44] "Peptides 4h_05"                      
## [45] "Peptides 4h_06"                      
## [46] "Peptides 24h_01"                     
## [47] "Peptides 24h_02"                     
## [48] "Peptides 24h_03"                     
## [49] "Peptides 24h_04"                     
## [50] "Peptides 24h_05"                     
## [51] "Peptides 24h_06"                     
## [52] "Protein IDs"                         
## [53] "Majority protein IDs"                
## [54] "Protein names"                       
## [55] "Gene names"                          
## [56] "id"

For prefiltering, one would first want to identify contaminants and bad IDs. There are three columns we should focus on: ‘Reverse’, ‘Potential Contaminant’, and ‘Only identified by site’.

Let’s take a look at what values are in these columns.

unique(data$Reverse)

## [1] NA

unique(data$`Potential contaminant`)

## [1] NA

unique(data$`Only identified by site`)

## [1] NA

Since this data was pre-filtered, there are no ‘hits’ for these rows. No rows have a flag indicating potential contamination or any other issue. If there were, it would look something like this: [1] "+".

Even though the bad rows have already been filtered out, here is what we would do if we did need to remove rows that were flagged. Assuming the "+" character is the flag, we would filter like this. The combination of ! (not) and %in% means to remove all rows equal to "+". At the end we print the percent of rows remaining after the filtering.

pre_contamination <- length(rownames(data))

data <- subset(data, !Reverse %in% "+")
data <- subset(data, !`Potential contaminant` %in% "+")
data <- subset(data, !`Only identified by site` %in% "+")

post_contamination <- length(rownames(data))

(pre_contamination/post_contamination) * 100

## [1] 100

As expected, we can see nothing was actually removed. Now let’s grab the information that we are actually interested in. We are using grep() to grab the first 18 columns or the intensity columns. The first argument to grep() is a regular expression that returns TRUE for any column name that begins with "0h_", "4h_", or "24h_", and FALSE otherwise.

intensity_cols <- grep("^(0h|4h|24h)_", colnames(data), value = TRUE)

analysis_data <- data %>%
  dplyr::select(all_of(intensity_cols))

head(analysis_data)

##      0h_01    0h_02    0h_03    0h_04    0h_05    0h_06    4h_01    4h_02
## 1 26.38970 25.84642 26.11164 24.26529 25.72419 26.07712 25.81389 26.13796
## 2 25.30655 32.27539 23.87797 32.47760 32.48539 24.57631 31.98132 24.57217
## 3 27.66662 27.07643 27.87050 27.48401 27.42021 27.67761 27.21011 28.01956
## 4 32.13081 31.79087 32.02677 31.91609 32.08689 31.82104 31.71959 31.97141
## 5 25.76018 26.22087 26.40334 24.64637 25.86446 26.11884 25.79823 26.07853
## 6 26.85792 27.07948 26.69112 26.70041 27.09678 26.68819 26.49919 26.62226
##      4h_03    4h_04    4h_05    4h_06   24h_01   24h_02   24h_03   24h_04
## 1 25.97062 25.78970 25.68664 26.40839 26.06528 25.80087 24.91407 23.87314
## 2 31.97641 25.26679 32.05093 24.74257 32.34056 31.95695 31.88965 32.16726
## 3 27.28542 27.72956 26.98414 28.03790 27.41844 27.34359 27.28824 26.89597
## 4 31.82487 31.86625 31.49348 31.76892 31.87339 31.62279 31.80355 31.61881
## 5 25.48025 26.21660 26.54041 26.10203 25.04800 26.74701 26.48800 25.41668
## 6 26.64623 26.30626 26.61975 26.26974 27.18773 27.04218 27.12119 26.85579
##     24h_05   24h_06
## 1 25.82260 24.62935
## 2 31.67745 31.80858
## 3 27.37227 25.06058
## 4 31.74859 31.58428
## 5 26.35795 26.42686
## 6 26.96521 27.22775

Look at the range (minimum and maximum) values.

range(analysis_data, na.rm = TRUE)

## [1] 22.64199 36.73819

The range is consistent with log₂ transformed data. For fun, let’s undo the log transformation and plot the resulting distribution.

raw_analysis_data <- 2^(analysis_data) - 1

for (i in 1:ncol(raw_analysis_data)) {
  hist(raw_analysis_data[[i]], main=colnames(raw_analysis_data)[i])
}

As we can see, the distributions are all heavily right-skewed. Let’s observe what happens after the log2 transformation…

log2_analysis_data <- log2(raw_analysis_data + 1)

for (i in 1:ncol(log2_analysis_data)) {
  hist(log2_analysis_data[[i]], main=colnames(log2_analysis_data)[i])
}

While the distributions are still right-skewed, they are more normally distributed than before.

Additionally it is a good idea to remove all rows (proteins) with an abundance of zero across all samples, or even ones with low abundance across all samples even if they are not exactly zero, but as we will see it has been pre-filtered. There are no rows with all zero values.

filter <- raw_analysis_data[rowSums(raw_analysis_data) == 0, ]

nrow(filter)

## [1] 0

Process and standardize data

Confirm that all intensity columns are in numeric format.

mat <- analysis_data %>% dplyr::select(all_of(intensity_cols))

sapply(mat, class)

##     0h_01     0h_02     0h_03     0h_04     0h_05     0h_06     4h_01     4h_02 
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" 
##     4h_03     4h_04     4h_05     4h_06    24h_01    24h_02    24h_03    24h_04 
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" 
##    24h_05    24h_06 
## "numeric" "numeric"

Extract the gene name column from the original data and correct any duplicates. Assign those names as row names of the intensity data frame. Transpose the data frame, which converts it into a matrix with samples as rows and proteins as columns.

gene_names <- data$`Gene names`

rownames(mat) <- make.unique(gene_names)

mat_t <- t(mat)

Use the scale() function to standardize the matrix. By default, this is done on a column basis. That means each protein is standardized, or z-transformed, so that its mean is 0 and standard deviation 1. We do this so that we can compare the relative difference in intensity between proteins that may have very different absolute intensities.

mat_t_scaled <- scale(mat_t)

Exploratory PCA

Perform a principal component analysis (PCA) on the scaled matrix. PCA is the most basic multivariate analysis technique. It finds a set of axes in p-dimensional space, where p is the number of proteins, that explains the most variation in the data. PCA is often a good first step when you have multivariate data like omics data, because it is very quick to compute, does a good job of picking out broad patterns in the data, and can identify any weird or problematic samples.

pca <- prcomp(mat_t_scaled)

Calculate the variance explained across the first two principal components and build a dataframe that includes the sample IDs and the PCA scores for each of the samples on the first two components. Then plot the samples’ coordinates on the first two PCA axes, labeling the points.

var_exp <- (pca$sdev^2) / sum(pca$sdev^2)

pc1_pct <- round(var_exp[1] * 100, 1)
pc2_pct <- round(var_exp[2] * 100, 1)


pca_df <- data.frame(
  sample = rownames(pca$x),
  group  = sub("_.*", "", rownames(pca$x)),
  PC1 = pca$x[,1],
  PC2 = pca$x[,2]
)

ggplot(pca_df, aes(PC1, PC2, label = sample, color = group)) +
  geom_point(size = 4) +
  geom_text(vjust = -1, color = "black") +
  xlab(paste0("PC1 (", pc1_pct, "%)")) +
  ylab(paste0("PC2 (", pc2_pct, "%)")) +
  scale_color_manual(values = custom_colors) +
  theme_minimal()

The plot shows that the first two PCA axes do a good job separating the time points from each other. It looks very similar to the transcriptomics plot: the first axis separates the 4h time point from the other two, and the second axis separates 0h and 24h from each other.

Differential expression analysis

Set up design matrix

We will create a design matrix for the three timepoints (0h, 4h, 24h), as we did for the transcriptome in Lesson 1. This is just a matrix of 0s and 1s indicating which time point each of the 18 samples belongs to. We also define the pairwise comparisons we want in the object contrast_matrix; again this is the same as in Lesson 1.

samples <- colnames(mat)

groups <- sub("_.*", "", samples)
groups <- factor(groups, levels = c("0h", "4h", "24h"))


design <- model.matrix(~ 0 + groups)
colnames(design) <- c("H0", "H4", "H24")
design

##    H0 H4 H24
## 1   1  0   0
## 2   1  0   0
## 3   1  0   0
## 4   1  0   0
## 5   1  0   0
## 6   1  0   0
## 7   0  1   0
## 8   0  1   0
## 9   0  1   0
## 10  0  1   0
## 11  0  1   0
## 12  0  1   0
## 13  0  0   1
## 14  0  0   1
## 15  0  0   1
## 16  0  0   1
## 17  0  0   1
## 18  0  0   1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$groups
## [1] "contr.treatment"

contrast_matrix <- makeContrasts(
    H0_vs_H4  = H0  - H4,
    H0_vs_H24 = H0 - H24,
    H4_vs_H24 = H4 - H24,
    levels = design
)

Fit linear model

Use the function lmFit() from the limma package to fit a linear model to each protein (column of the expression matrix mat).

fit <- lmFit(mat, design)

Apply the defined contrasts to the linear‑model fit and then compute empirical Bayes moderated statistics for differential expression.

fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

Visualize results

Volcano plot

Extract the differential‑expression result tables for each of the three contrasts (H0_vs_H4, H0_vs_H24, and H4_vs_H24) using topTable().

res_H0_vs_H4  <- topTable(fit2, coef = "H0_vs_H4",  number = Inf)
res_H0_vs_H24 <- topTable(fit2, coef = "H0_vs_H24", number = Inf)
res_H4_vs_H24 <- topTable(fit2, coef = "H4_vs_H24", number = Inf)

res_H0_vs_H4$Gene  <- rownames(res_H0_vs_H4)
res_H0_vs_H24$Gene <- rownames(res_H0_vs_H24)
res_H4_vs_H24$Gene <- rownames(res_H4_vs_H24)

Create a volcano plot visualizing the differential‑expression results for the “0h vs 4h” comparison. The proteins with a significant positive fold change were more abundant at the 0h time point, and the ones with a significant negative fold change were more abundant at the 4h time point.

volcano_plot(res_H0_vs_H4, "0h vs 4h")

Do the same for the other two pairs of time points.

volcano_plot(res_H0_vs_H24, "0h vs 24h")

volcano_plot(res_H4_vs_H24, "4h vs 24h")

Venn diagram

The following code creates a Venn diagram showing how many proteins were differentially abundant between each pair of time points. First, we set a threshold that any pair of proteins with a log₂ fold change of at least |0.5| (a difference of \(2^{\frac{{1}{2}}}\) or ~ 1.4×) and an empirical Bayes adjusted p-value of 0.05 or less are considered differentially abundant. Then the next few lines of code generate a matrix with 3 columns, one for each pair of time points, filled with 0s for pairs of time points where a particular gene was not differentially abundant, and 1s where it was.

adjP_thr <- 0.01
logFC_thr <- 0.5

get_up_genes <- function(df) {
  rownames(df)[
    abs(df$logFC)  > logFC_thr &
    df$adj.P.Val < adjP_thr
  ]
}

up_H0_vs_H4  <- get_up_genes(res_H0_vs_H4)
up_H0_vs_H24 <- get_up_genes(res_H0_vs_H24)
up_H4_vs_H24 <- get_up_genes(res_H4_vs_H24)

up_mat <- sapply(list(up_H0_vs_H4, up_H0_vs_H24, up_H4_vs_H24),
                 function(x) ifelse(row.names(mat) %in% x, 1, 0))

colnames(up_mat) <- c("H0_vs_H4", "H0_vs_H24", "H4_vs_H24")

vennDiagram(up_mat, circle.col = custom_colors)

The Venn diagram above shows that four proteins had different abundance at all three time points. Here are their names.

shared_up_genes <- Reduce(intersect,
                          list(up_H0_vs_H4,
                               up_H0_vs_H24,
                               up_H4_vs_H24))

shared_up_genes

## [1] "Plin2"   "Impdh1"  "Itm2b"   "Timm17a"

UpSet plot

As we did in Lesson 1, you may also want to make an UpSet plot instead of or in addition to a Venn diagram. Here’s the code to do that. We exclude the proteins that were not differentially abundant between any pair of time points. We see that about 100 proteins were differentially abundant between only one of the three possible pairs of time points, and a smaller number are overlapping (the same information conveyed by the Venn diagram).

grays <- c("gray30", "gray15", "gray0")

nonzero_comb_matrix <- make_comb_mat(up_mat[rowSums(up_mat) > 0, ])
UpSet(nonzero_comb_matrix, pt_size = unit(5, "mm"),
      bg_col = "#F0F0F0", lwd = 3,
      column_title = "Overlapping differential proteins",
      comb_col = grays[comb_degree(nonzero_comb_matrix)],
      right_annotation = upset_right_annotation(nonzero_comb_matrix, 
                                                gp = gpar(col = custom_colors, fill = custom_colors),
                                                annotation_name_side = "top",
                                                axis_param = list(side = "top")))

Multiomics workshop part 2: Proteome data

Matt Burnham, Quentin D. Read, & Matthias Anagho-Mattanovich

July 07, 2026