---
title: "Untitled"
author: "Swanitha"
date: "2025-09-03"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Load Libraries
```{r libraries}
library("limma")
library("ggplot2")
library("readr")
library("mitch")
```

#Load Data
```{r load_data}
load("z_log10.RData")
```

## Sample Setup

```{r sample_sheet}
# Sample sheet setup
sample_sheet <- data.frame(
    Sample_ID = colnames(z_log10),
    Fraction = ifelse(grepl("Chrom", colnames(z_log10)), "Chromatin", "Soluble"),
    Replicate = c(1, 2, 3, 4, 1, 2, 3, 4),
      stringsAsFactors = FALSE
)

sample_sheet$Fraction <- as.factor(sample_sheet$Fraction)

print(sample_sheet)
# Check that Fraction is now a factor
str(sample_sheet$Fraction) 
```

## Differential Expression Analysis

```{r limma_analysis}
# Design matrix
design <- model.matrix(~ Fraction, data = sample_sheet)
print(design)  

# Limma analysis
fit <- lmFit(z_log10, design)
fit <- eBayes(fit)

# Results summary
summary(decideTests(fit))
```
Analysis: Overall Differential Expression
Significantly altered proteins: 5,644 (67.4% of total)
Upregulated in Soluble: 3,367 proteins (40.2%)
Downregulated in Soluble: 2,277 proteins (27.2%)
Non-significant: 2,734 proteins (32.6%)

## Extract Results

```{r results}
# Extract all results
results <- topTable(fit, coef = 2, number = Inf)

# Basic statistics
summary(results$adj.P.Val < 0.05)
summary(results$logFC)
```
Analysis: normalized data with subtle but statistically significant differences
```{r}
# Top 20 results
head(results, 20)

#Bottom 20 results
tail(results, 20)
```
Analysis: Most Significantly Downregulated in Soluble (Chromatin-enriched)
PUM3 (RNA-binding protein) - adj.P < 3.4 × 10⁻⁷
DDX10 (RNA helicase) - adj.P < 3.4 × 10⁻⁷
PDCD11 (ribosome biogenesis) - adj.P < 3.4 × 10⁻⁷
DKC1 (telomerase complex) - adj.P < 3.4 × 10⁻⁷
RSL1D1 (ribosome biogenesis) - adj.P < 3.4 × 10⁻⁷

Minimal Changes (P-values > 0.97):

STK25 (serine/threonine kinase) - P = 0.9998, essentially no difference
ATG101 (autophagy protein) - P = 0.9975
HES1 (transcription factor) - P = 0.9974
DNAJC21 (chaperone) - P = 0.9945

Functional Categories:
Metabolic enzymes: MLYCD, IDUA, KDSR
Transcription factors: E2F3, HES1
Protein quality control: DNAJC8, DNAJC21
Membrane proteins: TMEM167A, TM7SF3
Kinases: STK25, PI4K2B

## Classify Results

```{r classification}
# Significant proteins
upregulated <- results[results$logFC > 0 & results$adj.P.Val < 0.05, ]
downregulated <- results[results$logFC < 0 & results$adj.P.Val < 0.05, ]

# Count significant proteins
sig_up <- nrow(upregulated)
sig_down <- nrow(downregulated)

print(paste("Upregulated proteins:", sig_up))
print(paste("Downregulated proteins:", sig_down))
```

## Histone Analysis

```{r histones}
# Extract histone proteins
histones <- results[grepl("_H1|_H2A|_H2B|_H3|_H4", rownames(results)), ]
print(histones)
```

## Basic Visualization

```{r volcano_plot}
plot(results$logFC, -log10(results$P.Value), 
     xlab = "log2(Fold Change)", ylab = "-log10(P-value)",
     main = "Volcano Plot: Soluble vs Chromatin",
     col = "black", pch = 16)

# Add significance threshold line
abline(h = -log10(0.05), col = "gray", lty = 2)
abline(v = 0, col = "gray", lty = 2)

# Identify and highlight significant points
sig <- subset(results, P.Value < 0.05)
points(sig$logFC, -log10(sig$P.Value), col = "blue", pch = 16)

```

## Gene Symbols

```{r gene_symbols}
# Extract gene symbols
upregulated_genes <- gsub(".*_", "", rownames(upregulated))
downregulated_genes <- gsub(".*_", "", rownames(downregulated))

head(upregulated_genes, 10)
head(downregulated_genes, 10)
```

## Save Results

```{r save_results}
# Save to CSV files
write.csv(results, "differential_expression_results.csv", row.names = TRUE)
write.csv(results[results$adj.P.Val < 0.05, ], "significant_proteins.csv")
write.csv(histones, "histone_analysis.csv")
```

## Pathway analysis

```{r,pathway}

if ( ! file.exists("ReactomePathways_2025-09-02.gmt") ) {
  message("Reactome file not found")
  download.file("https://ziemann-lab.net/public/proteomics/prot_de_eda/ReactomePathways_2025-09-02.gmt")
} else {
  message("Reactome file found")
}

reactome <- gmt_import("ReactomePathways_2025-09-02.gmt")

res2 <- results[which(lapply(strsplit(rownames(results),"_"),length)==2),]

which(duplicated(sapply(strsplit(rownames(res2),"_"),"[[",2)))

rownames(res2) <- sapply(strsplit(rownames(res2),"_"),"[[",2)

m <- mitch_import(x=res2,DEtype="limma")

mres <- mitch_calc(x=m,genesets=reactome,minsetsize=5,cores=4,priority="effect")

dim(mres$enrichment_result)

head(subset(mres$enrichment_result,s.dist>0),20) %>%
  kbl(caption="Reactome pathway analysis (soluble fraction sets by effect size)") %>%
  kable_paper("hover", full_width = F)

head(subset(mres$enrichment_result,s.dist<0),20) %>%
  kbl(caption="Reactome pathway analysis (chromatin fraction sets by effect size)") %>%
  kable_paper("hover", full_width = F)

```
