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

#### 1. Load Required Packages
# Install required packages
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

# Install packages from CRAN
cran_packages <- c("pheatmap", "ggplot2", "dplyr", "ggrepel", "visNetwork")
install.packages(cran_packages)

# Install packages from Bioconductor
bioc_packages <- c("AnnotationDbi", "org.Hs.eg.db", "clusterProfiler", "STRINGdb")
BiocManager::install(bioc_packages)

### NOTE: For those who want to create a separate environment in RStudio for analysis, run the codes as below (same for DDA analysis): 
# Install the 'reticulate' package if not already installed
install.packages("reticulate")

# Load the reticulate package
library(reticulate)

# Install Miniconda via reticulate
install_miniconda()

# Check if conda is available
conda_binary()

# Create the environment using the provided YAML file
conda env create -f Supplemental_File_3.yaml

### load packages
library(pheatmap)       
library(ggplot2)        
library(dplyr)          
library(AnnotationDbi)  
library(org.Hs.eg.db)   
library(clusterProfiler) 
library(STRINGdb)       
library(visNetwork)     
library(ggrepel)        

#### 2. Data Loading and Preprocessing with Proteotypic Peptide Filter
# Read protein group matrix
result <- read.table("~/PDAC_report.pg_matrix.tsv",
                     stringsAsFactors = FALSE, 
                     sep = '\t', 
                     row.names = 1, 
                     header = TRUE)

# Apply proteotypic peptide filter (N.Proteotypic.Sequences ≥ 2)
filtered_data <- result[result$N.Proteotypic.Sequences >= 2, ]

# Extract expression matrix (columns 6-11) and remove missing values
prot_mat_withNA <- filtered_data[, 6:11]
prot_mat <- na.omit(prot_mat_withNA)

#### 3. Median Normalization
median_normalize <- function(mat) {
  col_medians <- apply(mat, 2, median, na.rm = TRUE)
  normalized_mat <- t(t(mat) / col_medians)
  return(normalized_mat)
}
prot_mat_normalized <- median_normalize(prot_mat)

#### 4. Group-wise Analysis
# Split into Normal (first 3) and PDAC (last 3) groups
protein_mat_normal <- prot_mat_normalized[, 1:3]
protein_mat_PDAC <- prot_mat_normalized[, 4:6]

# Calculate mean expression
PDAC_mean <- rowMeans(protein_mat_PDAC, na.rm = TRUE)
normal_mean <- rowMeans(protein_mat_normal, na.rm = TRUE)

# Perform t-tests between groups
PDAC_pvalue <- apply(prot_mat_normalized, 1, function(x) {
  t.test(x[4:6], x[1:3])$p.value
})

#### 5. Quality Control Visualization
# Log2FC distribution
hist(log2(PDAC_mean / normal_mean),
     main = "Log2 Fold Change Distribution (Proteotypic ≥ 2)",
     xlab = "Log2(PDAC/Normal)")

# Volcano plot (basic)
plot(log2(PDAC_mean/normal_mean), -log10(PDAC_pvalue),
     main = "Volcano Plot (Proteotypic ≥ 2)",
     xlab = "Log2 Fold Change",
     ylab = "-Log10(p-value)",
     pch = 19, col = ifelse(abs(log2(PDAC_mean/normal_mean)) > 1 & PDAC_pvalue < 0.05, "red", "black"))
abline(h = -log10(0.05), col = "blue", lty = 2)
abline(v = c(-1, 1), col = "blue", lty = 2)

#### 6. Enhanced Volcano Plot (ggplot2)

# Create the data frame
results_df <- data.frame(
  ProteinID = rownames(prot_mat_normalized),
  Gene = mapIds(org.Hs.eg.db, 
                keys = rownames(prot_mat_normalized),
                column = "SYMBOL",
                keytype = "UNIPROT",
                multiVals = "first"),
  Log2FC = log2(PDAC_mean / normal_mean),
  p.value = PDAC_pvalue,
  N.Proteotypic = filtered_data[rownames(prot_mat_normalized), "N.Proteotypic.Sequences"]
) %>% 
  mutate(
    # Correct p-values for multiple testing
    FDR = p.adjust(p.value, method = "BH"),
    
    # Significance classification
    Significance = case_when(
      FDR < 0.05 & abs(Log2FC) > 0.6 ~ "FDR < 0.05 & |FC| > 0.6",
      p.value < 0.05 & abs(Log2FC) > 0.6 ~ "Nominal p < 0.05 & |FC| > 0.6",
      TRUE ~ "Not Significant"
    )
  )

# Check for missing gene names
sum(is.na(results_df$Gene))

# Replace missing gene names with ProteinID
results_df$Gene[is.na(results_df$Gene)] <- results_df$ProteinID[is.na(results_df$Gene)]

# Create the volcano plot
ggplot(results_df, aes(x = Log2FC, y = -log10(p.value))) +
  geom_point(
    aes(color = Significance, size = N.Proteotypic), 
    alpha = 0.6
  ) +
  scale_color_manual(
    values = c(
      "FDR < 0.05 & |FC| > 0.6" = "red",
      "Nominal p < 0.05 & |FC| > 0.6" = "orange", 
      "Not Significant" = "grey80"
    ),
    labels = c(
      "FDR < 0.05 & |FC| > 0.6" = "FDR < 0.05 & |log2FC| > 0.6",
      "Nominal p < 0.05 & |FC| > 0.6" = "Nominal p < 0.05 & |log2FC| > 0.6",
      "Not Significant" = "Not Significant"
    )
  ) +
  geom_text_repel(
    data = subset(results_df, Significance == "FDR < 0.05 & |FC| > 0.6" | Significance == "Nominal p < 0.05 & |FC| > 0.6"),
    aes(label = Gene),
    size = 3,
    max.overlaps = 20,
    min.segment.length = 0.2,
    box.padding = 0.5,
    segment.color = "grey50"
  ) +
  geom_hline(
    yintercept = -log10(0.05), 
    linetype = "dashed", 
    color = "blue"
  ) +
  geom_vline(
    xintercept = c(-0.6, 0.6), 
    linetype = "dashed", 
    color = "blue"
  ) +
  labs(
    title = "PDAC Serum Tumor vs Normal (Proteotypic ≥ 2 peptides)",
    subtitle = "Median-normalized data; Red: FDR < 0.05 & |log2FC| > 0.6; Orange: Nominal p < 0.05 & |log2FC| > 0.6",
    x = "log2(Fold Change)",
    y = "-log10(p-value)",
    color = "Significance",
    size = "# Proteotypic Peptides"
  ) +
  theme_classic() +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold")
  )
#### 7. Heatmap of Significant Proteins
# Get significant proteins
sig_proteins <- results_df %>% 
  filter(Significance == "FDR < 0.05 & |FC| > 0.6" | Significance == "Nominal p < 0.05 & |FC| > 0.6") %>% 
  pull(ProteinID)
sig_expr <- prot_mat_normalized[sig_proteins, ]

# Z-score normalization
sig_expr_z <- t(scale(t(sig_expr)))

# Custom column names
colnames(sig_expr_z) <- c("Normal_Serum_2112", "Normal_Serum_2195", "Normal_Serum_2377",
                          "PDAC_Serum_8526", "PDAC_Serum_8568", "PDAC_Serum_8696")

# Plot heatmap
pheatmap(sig_expr_z,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         border_color = NA,
         cluster_rows = TRUE,
         cluster_cols = FALSE,
         show_rownames = FALSE,
         annotation_col = data.frame(
           Group = factor(c(rep("Normal", 3), rep("PDAC", 3))),
           row.names = colnames(sig_expr_z)
         ),
         main = "Significant Proteins Heatmap (Proteotypic ≥ 2)")

#### 8. Functional Enrichment Analysis
# Get Entrez IDs
entrez_ids <- mapIds(org.Hs.eg.db, keys = sig_proteins,
                     column = "ENTREZID", keytype = "UNIPROT", multiVals = "first")

# GO Enrichment (Biological Process)
go_enrich <- enrichGO(gene = na.omit(entrez_ids),
                      OrgDb = org.Hs.eg.db,
                      ont = "BP",
                      pvalueCutoff = 0.05,
                      readable = TRUE)

dotplot(go_enrich, showCategory = 20, 
        title = "GO Enrichment (Proteotypic ≥ 2)",
        label_format = 30) +
  theme(axis.text.y = element_text(size = 10, margin = margin(r = 15)))

# KEGG Pathway Enrichment
kegg_enrich <- enrichKEGG(gene = na.omit(entrez_ids),
                          organism = "hsa",
                          pvalueCutoff = 0.05)

dotplot(kegg_enrich, showCategory = 20, 
        title = "KEGG Pathway Enrichment (Proteotypic ≥ 2)")

#### 9. Protein-Protein Interaction Network
# STRINGdb setup
string_db <- STRINGdb$new(version = "11.5", species = 9606, score_threshold = 400)

# Map proteins to STRING IDs
string_mapped <- string_db$map(data.frame(protein = sig_proteins), "protein")

# Generate network
interactions <- string_db$get_interactions(string_mapped$STRING_id)

# Prepare nodes and edges
nodes <- data.frame(
  id = string_mapped$STRING_id,
  label = string_mapped$protein,
  title = mapIds(org.Hs.eg.db, keys = string_mapped$protein,
                 column = "SYMBOL", keytype = "UNIPROT", multiVals = "first"),
  value = filtered_data[string_mapped$protein, "N.Proteotypic.Sequences"] # Size by proteotypic count
)

edges <- interactions %>% mutate(
  from = from, to = to, width = combined_score/200)

# Interactive visualization
visNetwork(nodes, edges) %>%
  visNodes(size = "value", # Scale node size by proteotypic count
           font = list(size = 12)) %>%
  visEdges(smooth = TRUE) %>%
  visPhysics(solver = "forceAtlas2Based",
             forceAtlas2Based = list(gravitationalConstant = -50)) %>%
  visOptions(highlightNearest = list(enabled = TRUE, hover = TRUE),
             nodesIdSelection = TRUE) %>%
  visLayout(randomSeed = 123) # For reproducible layout