### DDA Proteomics Data Analysis Pipeline (with Proteotypic Peptide Filter)

# Install required packages (run once)
install.packages(c("tidyverse", "limma", "EnhancedVolcano", "clusterProfiler", "pheatmap", "RColorBrewer"))
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db")

# Load required packages
library(tidyverse)
library(limma)
library(EnhancedVolcano)
library(clusterProfiler)
library(org.Hs.eg.db)
library(pheatmap)
library(RColorBrewer)

# 1. Data Loading and Preparation ------------------------------------------
protein_data <- read_tsv("~/combined_protein.tsv") %>% 
  as.data.frame()

# 2. Filter for proteins with Combined Unique Spectral Count ≥ 2 -----------
protein_data <- protein_data %>% 
  filter(`Combined Unique Spectral Count` >= 2)

# 3. Extract Expression Matrix --------------------------------------------
maxlfq_cols <- grep("MaxLFQ", colnames(protein_data))
exp_matrix <- as.matrix(protein_data[, maxlfq_cols])
rownames(exp_matrix) <- protein_data$`Protein ID`

# Remove rows with NA values
exp_matrix <- exp_matrix[complete.cases(exp_matrix), ]

# 4. Create Sample Groups -------------------------------------------------
# Automatically detect groups from column names (P = Normal, T = Tumor)
colnames(exp_matrix) <- gsub("20150127_liver_XH03_|_MaxLFQ Intensity", "", colnames(exp_matrix))
group <- factor(ifelse(grepl("^P_", colnames(exp_matrix)), "Normal", "Tumor"))

# 5. Quantile normalization
normalized_exp <- normalizeBetweenArrays(exp_matrix, method = "quantile")

# Visualize normalization effect
par(mfrow = c(1, 2))
boxplot(exp_matrix, main = "Before Normalization", las = 2)
boxplot(normalized_exp, main = "After Quantile Normalization", las = 2)

# 6. Calculate Group Means ------------------------------------------------
normal_mean <- rowMeans(normalized_exp[, group == "Normal"])
tumor_mean <- rowMeans(normalized_exp[, group == "Tumor"])

# 7. Calculate log2 Fold Change -------------------------------------------
log2FC <- log2((tumor_mean) / (normal_mean))  # Added pseudocount to avoid division by zero

# 8. Plot Fold Change Distribution ----------------------------------------
hist(log2FC,
     main = "Log2 Fold Change Distribution (Quantile Normalized)",
     xlab = "Log2(Tumor/Normal)",
     col = "lightblue",
     breaks = 30)

# 9. Differential Expression Analysis (t-test) ---------------------------
tumor_pvalue <- sapply(1:nrow(normalized_exp), function(i) {
  x <- normalized_exp[i, group == "Tumor"]
  y <- normalized_exp[i, group == "Normal"]
  
  # Check for zero variance or insufficient samples
  if (sd(x) == 0 || sd(y) == 0 || length(x) < 2 || length(y) < 2) {
    return(NA)
  } else {
    t.test(x, y)$p.value
  }
})

# Create results dataframe with multiple testing correction
diff_results <- data.frame(
  Protein.ID = rownames(normalized_exp),
  Mean_Normal = normal_mean,
  Mean_Tumor = tumor_mean,
  logFC = log2FC,
  p.value = tumor_pvalue,
  adj.P.Val = p.adjust(tumor_pvalue, method = "BH"), # Benjamini-Hochberg correction
  stringsAsFactors = FALSE
) %>% 
  left_join(
    protein_data[, c("Protein ID", "Gene", "Description", "Combined Unique Spectral Count")], 
    by = c("Protein.ID" = "Protein ID") 
  )

# Filter out invalid rows (NA or infinite values)
diff_results_filtered <- diff_results %>%
  filter(!is.na(logFC) & !is.infinite(logFC) &
           !is.na(p.value) & !is.infinite(p.value))

# 10. Generate Volcano Plot -----------------------------------------------
EnhancedVolcano(diff_results_filtered,
                lab = diff_results_filtered$Gene,
                x = 'logFC',
                y = 'adj.P.Val', # Using adjusted p-values or y = 'p.value'
                pCutoff = 0.05,
                FCcutoff = 1,
                title = "HCC Tumor vs Normal (t-test)",
                subtitle = "Quantile-normalized data | FDR-adjusted p-values",
                legendPosition = 'right',
                drawConnectors = TRUE,
                widthConnectors = 0.5,
                max.overlaps = 100,
                labSize = 3,
                pointSize = 1.5,
                ylim = c(0, max(-log10(diff_results_filtered$adj.P.Val), na.rm = TRUE) * 1.1))

# 11. Generate Heatmap ----------------------------------------------------
top_diff <- diff_results_filtered %>% 
  arrange(adj.P.Val) %>% 
  head(50) %>% 
  pull(Protein.ID)

heatmap_data <- normalized_exp[top_diff, ]

# Create sample annotation
annotation_col <- data.frame(
  Condition = group,
  row.names = colnames(heatmap_data)
)

# Custom color palette
heatmap_colors <- colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100)

pheatmap(heatmap_data,
         scale = "row",
         annotation_col = annotation_col,
         show_rownames = FALSE,
         main = "Top 50 Differentially Expressed Proteins in HCC",
         color = heatmap_colors,
         fontsize_col = 8,
         clustering_method = "ward.D2")

# 12. Pathway Enrichment Analysis ----------------------------------------
# Get significant genes (adj.P.Val < 0.05 and |logFC| > 1)
sig_genes <- diff_results_filtered %>% 
  filter(adj.P.Val < 0.05 & abs(logFC) > 1) %>% 
  pull(Gene) %>% 
  str_split(";") %>% 
  unlist() %>% 
  unique()

# Convert gene symbols to ENTREZ IDs
entrez_ids <- bitr(sig_genes, 
                   fromType = "SYMBOL",
                   toType = "ENTREZID",
                   OrgDb = org.Hs.eg.db)

# GO enrichment analysis
go_enrich <- enrichGO(gene = entrez_ids$ENTREZID, 
                      OrgDb = org.Hs.eg.db, 
                      ont = "BP",
                      pvalueCutoff = 0.05, # More stringent cutoff
                      qvalueCutoff = 0.2,
                      readable = TRUE)

# KEGG pathway analysis
kegg_enrich <- enrichKEGG(gene = entrez_ids$ENTREZID, 
                          organism = 'hsa',
                          pvalueCutoff = 0.05,
                          qvalueCutoff = 0.2)

# Visualize enrichment results
if (nrow(go_enrich@result) > 0) {
  dotplot(go_enrich, 
          showCategory = 15,
          title = "GO Biological Process Enrichment (FDR < 0.05)",
          font.size = 10,
          label_format = 30) +
    theme(axis.text.y = element_text(size = 10))
} else {
  message("No significant GO terms found.")
}

if (nrow(kegg_enrich@result) > 0) {
  dotplot(kegg_enrich, 
          showCategory = 15,
          title = "KEGG Pathway Enrichment (FDR < 0.05)",
          font.size = 10)
} else {
  message("No significant KEGG pathways found.")
}

# 13. Save Results -------------------------------------------------------
write_csv(diff_results_filtered, "HCC_Differential_Expression_Results.csv")