############################################################################################
#### MULTIOMICS WORKSHOP 2026                                                           ####
#### LESSON 2: PROTEOME DATA                                                            ####
#### =======================                                                            ####
#### This worksheet contains an incomplete version of the code presented in the lesson. ####
#### Fill in all ... with code.                                                         ####
############################################################################################
############################################################################################


### Setup

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

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

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"))

}

# Path to raw data file. 
# If you aren't running this on the demo server, uncomment the url() line and run it to read from GitHub Pages URL.
raw <- "data/Imputed_proteinGroups_perseus.txt"
# raw <- url("https://usda-ree-ars.github.io/SEAStatsData/omics2026/Imputed_proteinGroups_perseus.txt")

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

### Inspect data

colnames(data)

unique(data$Reverse)
unique(data$`Potential contaminant`)
unique(data$`Only identified by site`)

### Remove flagged rows

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

### Get data columns

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

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

head(analysis_data)

range(analysis_data, na.rm = TRUE)

### Plot untransformed and log-transformed data

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])
}

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])
}

### Remove zero-expression rows

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

nrow(filter)

### Prepare data matrix

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

sapply(mat, class)

gene_names <- data$`Gene names`

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

mat_t <- t(mat)

### Standardize data, do PCA, and plot

mat_t_scaled <- scale(mat_t) 

pca <- prcomp(mat_t_scaled)

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()

### Create design and contrast matrices

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

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

### Differential abundance analysis, with empirical Bayes correction

fit <- lmFit(mat, design)

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

### Extract results and create Volcano plots

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)

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

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

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

### Venn Diagram

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)

### View differentially abundant proteins

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

shared_up_genes

### UpSet plot

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")))