---
title: "Proteomics EDA"
author: "Swanitha"
date: "2025-08-24"
output: html_document
---

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

```{r}
# Load required libraries
library("pcaMethods")
library("limma")

```

```{r}
# Load the data
# PATH projects/proteomics/workingdir02/dia-quant-output/report.pg_matrix.tsv
x <- read.table("report.pg_matrix.tsv", header=TRUE, sep="\t")

# Clean column names
colnames(x) <- gsub("X.home.projects.proteomics.PXD058340.20240911_AST0_NEO4_IAH_collab_SAG_U2OS_Frac_Control_","", colnames(x))
colnames(x) <- gsub("_uncalibrated.mzML", "", colnames(x))

# Set row names
rownames(x) <- paste(x$Protein.Ids, x$Genes, sep = "_")

# Extract intensity data
y <- x[, 6:13]

# Order columns alphabetically
y <- y[, order(colnames(y))]

# Display data dimensions and first few rows
cat("Data dimensions:", nrow(y), "proteins x", ncol(y), "samples\n")
cat("FIRST 6 ROWS OF INTENSITY DATA\n")
print(head(y))
```

```{r}
# Missing Data Analysis
# Overall missing data statistics
total_values <- nrow(y) * ncol(y)
total_missing <- sum(is.na(y))
missing_percentage <- (total_missing / total_values) * 100

cat("\nOVERALL MISSING DATA SUMMARY\n")
cat("Total values:", total_values, "\n")
cat("Total missing:", total_missing, "\n")
cat("Overall missing percentage:", round(missing_percentage, 2), "%\n")

# Missing data per sample
missing_per_sample <- colSums(is.na(y))
missing_per_sample_pct <- (missing_per_sample / nrow(y)) * 100

missing_summary <- data.frame(
  Sample = names(missing_per_sample),
  Missing_Count = missing_per_sample,
  Missing_Percent = round(missing_per_sample_pct, 2)
)

cat("\nMISSING DATA PER SAMPLE\n")
print(missing_summary)
```

```{r}
# Missing data per protein
nas <- apply(X = y, MARGIN = 1, function(x){ 
  length(which(is.na(x))) 
})
missing_proteins_table <- table(nas)

cat("\nDATA DISTRIBUTION ACROSS PROTEINS\n")
cat("Number of proteins by missing value count:\n")
print(missing_proteins_table)

cat("\nSummary of missing values per protein:\n")
print(summary(nas))
```
Analysis
7.55% overall missing data, sample-wise missing data ranges from 6.64%(sol_01) to 8.96% (Chrom_01). The pattern suggests data could be missing completely at random (MCAR)
```{r}
# Raw Data Distribution Analysis
# Basic statistics for each sample
raw_stats <- data.frame(
  Sample = colnames(y),
  Min = apply(y, 2, min, na.rm = TRUE),
  Q1 = apply(y, 2, quantile, 0.25, na.rm = TRUE),
  Median = apply(y, 2, median, na.rm = TRUE),
  Mean = apply(y, 2, mean, na.rm = TRUE),
  Q3 = apply(y, 2, quantile, 0.75, na.rm = TRUE),
  Max = apply(y, 2, max, na.rm = TRUE),
  SD = apply(y, 2, sd, na.rm = TRUE)
)

# Round only numeric columns
raw_stats[, 2:8] <- round(raw_stats[, 2:8], 2)

cat("\nRAW DATA STATISTICS PER SAMPLE\n")
print(raw_stats)

# Summary of all samples combined
cat("\nDATA SUMMARY\n")
print(summary(y))
```
Analysis
Chromatin samples show higher maximum values and variance
```{r}
# Data Visualization - Raw Distribution
# Boxplot of raw abundances
par(mfrow = c(1, 1), mar = c(5, 4, 3, 1))

# Raw scale boxplot
boxplot(y, las = 2, main = "Raw Abundance Distribution", 
        ylab = "Abundance", cex.axis = 0.6, col = rainbow(ncol(y)))

# Log-transformed boxplot for better visualization
y_log <- log10(y + 1)  # Add 1 to handle zeros if any
boxplot(y_log, las = 2, main = "Log10(Abundance + 1) Distribution", 
        ylab = "Log10(Abundance + 1)", cex.axis = 0.6, col = rainbow(ncol(y)))

par(mfrow = c(1, 1))
```

```{r}
# Protein Abundance Analysis
# Calculate protein statistics
protein_stats <- data.frame(
  Protein = rownames(y),
  Mean_Abundance = rowMeans(y, na.rm = TRUE),
  Median_Abundance = apply(y, 1, median, na.rm = TRUE),
  SD_Abundance = apply(y, 1, sd, na.rm = TRUE),
  CV = apply(y, 1, sd, na.rm = TRUE) / rowMeans(y, na.rm = TRUE),
  Missing_Count = nas
)

cat("\nPROTEIN ABUNDANCE SUMMARY\n")
cat("Mean abundance across all proteins:\n")
print(summary(protein_stats$Mean_Abundance))

cat("\nCoefficient of Variation (CV) summary:\n")
print(summary(protein_stats$CV))

# Top 10 most abundant proteins
high_abundance <- head(protein_stats[order(protein_stats$Mean_Abundance, decreasing = TRUE), ], 10)

cat("\nTOP 10 MOST ABUNDANT PROTEINS\n")
print(high_abundance[, c("Protein", "Mean_Abundance", "Missing_Count")])
```

```{r}
# Top 10 most variable proteins by CV
high_variable <- head(protein_stats[order(protein_stats$CV, decreasing = TRUE), ], 10)

cat("\nTOP 10 MOST VARIABLE PROTEINS (by Coefficient of Variation)\n")
#Round only the numeric columns, keep Protein column as-is
high_variable_display <- high_variable[, c("Protein", "Mean_Abundance", "CV", "Missing_Count")]
high_variable_display[, c("Mean_Abundance", "CV")] <- round(high_variable_display[, c("Mean_Abundance", "CV")], 3)
high_variable_display[, "Missing_Count"] <- round(high_variable_display[, "Missing_Count"], 0)  # Missing count should be integer
print(high_variable_display)
```

```{r}
# Visualization of CV distribution
hist(protein_stats$CV[!is.infinite(protein_stats$CV)], 
     breaks = 50, main = "Distribution of Coefficient of Variation", 
     xlab = "CV", ylab = "Frequency", 
     col = "lightblue", border = "darkblue")
```

```{r}
# Sample Correlation Matrix
# Calculate correlation matrix
cor_matrix <- cor(y, use = "complete.obs")

cat("\nSAMPLE CORRELATION MATRIX\n")
print(round(cor_matrix, 3))

# Identify sample groups
chrom_samples <- grepl("Chrom", colnames(y))
sol_samples <- grepl("Sol", colnames(y))

# Calculate group-wise correlations
chrom_cor <- cor_matrix[chrom_samples, chrom_samples]
sol_cor <- cor_matrix[sol_samples, sol_samples]
between_cor <- cor_matrix[chrom_samples, sol_samples]

cat("\nCORRELATION SUMMARY\n")
cat("Average within-Chromatin correlation:", round(mean(chrom_cor[upper.tri(chrom_cor)]), 3), "\n")
cat("Average within-Soluble correlation:", round(mean(sol_cor[upper.tri(sol_cor)]), 3), "\n")
cat("Average between-condition correlation:", round(mean(between_cor), 3), "\n")
```

```{r}
# Simple correlation heatmap using base R
heatmap(cor_matrix, main = "Sample Correlation Heatmap", 
        col = heat.colors(100), margins = c(8, 8))
```

```{r}
# Data Quality Assessment - Define z before PCA
# Filter proteins with fewer missing values (e.g., less than 3 missing)
z <- y[which(nas < 3), ]
cat("\nDATA QUALITY ASSESSMENT\n")
cat("Original data:", nrow(y), "proteins\n")
cat("After filtering (< 3 missing values):", nrow(z), "proteins\n")
cat("Proteins retained:", round(nrow(z)/nrow(y) * 100, 1), "%\n")
```

```{r}
#NORMALIZATION PIPELINE
cat("\nNORMALIZATION PIPELINE\n")

#Load required libraries
library(limma)

#Handle remaining missing values in filtered data
z_complete <- z
for(i in 1:ncol(z_complete)) {
  if(any(is.na(z_complete[,i]))) {
    min_val <- min(z_complete[,i], na.rm = TRUE)
    z_complete[is.na(z_complete[,i]), i] <- min_val * 0.1
  }
}

cat("Missing values after imputation:", sum(is.na(z_complete)), "\n")

#Log10 transformation
z_log10 <- log10(z_complete + 1)
cat("Log10 transformation completed\n")

#Between-array normalization
z_normalized <- normalizeBetweenArrays(z_log10, method = "quantile", targets = NULL, cyclic.method = "fast")
cat("Between-array normalization completed\n")

cat("Final normalized data dimensions:", nrow(z_normalized), "x", ncol(z_normalized), "\n")
```

```{r}
# NORMALIZATION VALIDATION
cat("\nNORMALIZATION VALIDATION\n")

# Before and after comparison
cat("Raw data range:", round(range(z_complete), 0), "\n")
cat("Normalized data range:", round(range(z_normalized), 2), "\n")

# Sample median variability
raw_medians <- apply(z_complete, 2, median)
norm_medians <- apply(z_normalized, 2, median)

cat("CV of sample medians (raw):", round(sd(raw_medians)/mean(raw_medians), 3), "\n")
cat("CV of sample medians (normalized):", round(sd(norm_medians)/mean(norm_medians), 3), "\n")

# Correlation preservation
cor_raw <- cor(z_complete)
cor_norm <- cor(z_normalized)

chrom_idx <- grepl("Chrom", colnames(z_normalized))
sol_idx <- grepl("Sol", colnames(z_normalized))

cat("Within-Chromatin correlation (normalized):", 
    round(mean(cor_norm[chrom_idx, chrom_idx][upper.tri(cor_norm[chrom_idx, chrom_idx])]), 3), "\n")
cat("Within-Soluble correlation (normalized):", 
    round(mean(cor_norm[sol_idx, sol_idx][upper.tri(cor_norm[sol_idx, sol_idx])]), 3), "\n")
```

```{r}
# VISUALIZATION COMPARISON
par(mfrow = c(2, 2), mar = c(6, 4, 3, 1))

# Raw data
boxplot(log10(z_complete + 1), las = 2, main = "Raw Data (Log10 scale)", 
        cex.axis = 0.6, col = rainbow(ncol(z_complete)))

# Log10 transformed
boxplot(z_log10, las = 2, main = "Log10 Transformed", 
        cex.axis = 0.6, col = rainbow(ncol(z_log10)))

# Final normalized(quantile)
boxplot(z_normalized, las = 2, main = "Final Normalized", 
        cex.axis = 0.6, col = rainbow(ncol(z_normalized)))

par(mfrow = c(1, 1))
```

```{r}
# Sample completeness
sample_completeness <- (1 - missing_per_sample_pct/100) * 100
cat("\nSample completeness (%):\n")
for(i in 1:length(sample_completeness)) {
  cat(names(sample_completeness)[i], ":", round(sample_completeness[i], 1), "%\n")
}
```

```{r}
# PCA Analysis using pcaMethods
library(pcaMethods)
cat("\nPCA ANALYSIS using pcaMethods\n")

# Use filtered data (z) with fewer missing values
# Transpose data so samples are rows (pcaMethods expects this format)
z_t <- t(z)

# Perform PPCA (handles missing values)
pc <- pca(z_t, nPcs = 3, method = "ppca")

# Display PCA results
print(pc)

cat("\nPCA Summary:\n")
cat("Number of components:", pc@nPcs, "\n")
cat("Method used:", pc@method, "\n")
cat("Variance explained by each PC:\n")
print(round(pc@R2, 3))

cat("\nCumulative variance explained:\n")
print(round(pc@R2cum, 3))

# Scores plot (samples)
par(mar = c(5,4,3,1))
slplot(pc, main = "PCA Scores Plot", cex = 0.8, cex.axis = 0.8)

# Manual scores plot with sample group colors
scores <- pc@scores
chrom_color <- ifelse(grepl("Chrom", rownames(scores)), "red", "blue")

plot(scores[,1], scores[,2], col = chrom_color, pch = 19, cex = 1.5,
     main = "PCA Scores - PC1 vs PC2",
     xlab = paste0("PC1 (", round(pc@R2[1]*100, 1), "% variance)"),
     ylab = paste0("PC2 (", round(pc@R2[2]*100, 1), "% variance)"))

# Add sample labels
text(scores[,1], scores[,2], labels = rownames(scores), pos = 3, cex = 0.7)

# Add legend
legend("topright", legend = c("Chromatin", "Soluble"), 
       col = c("red", "blue"), pch = 19, cex = 0.8)
```
