Dataset 1

# Load count matrix and metadata
counts_D1 <- read.delim("counts_D1", header = TRUE, row.names = 1)
metadata_D1 <- read.csv("metadata_D1.csv", header = TRUE)

# Check distribution of counts
total_counts_per_gene_D1 <- rowSums(counts_D1)
hist(total_counts_per_gene_D1, 
     breaks=100, 
     main="Distribution of Total Counts per Gene", 
     xlab="Total counts across all samples",
     ylab="Number of genes")
summary(total_counts_per_gene_D1)

# Filter low expressed genes
counts_D1 <- counts_D1[which(rowSums(counts_D1) > 50), ]

# Round counts for DESeq2
counts_D1 <- round(counts_D1)

# Create DESeq2 object
library(DESeq2)
dds_D1 <- DESeqDataSetFromMatrix(countData = counts_D1,
                                colData = metadata_D1,        
                                design = ~Treatment)

dds_D1

# Run DESeq2 analysis
dds_D1 <- DESeq(dds_D1)

# Variance stabilizing transformation
vsdata_D1 <- vst(dds_D1, blind = TRUE)

summary(dds_D1)
#Get VST matrix (genes as rows, samples as columns)
vst_matrix_D1 <- assay(vst(dds_D1, blind = TRUE))
head(vst_matrix_D1)

library(reshape2)
vst_df_D1 <- as.data.frame(vst_matrix_D1)
vst_df_D1$gene_id <- rownames(vst_df_D1)

vst_final_D1 <- melt(vst_df_D1, 
                 id.vars = "gene_id", 
                 variable.name = "sample_id", 
                 value.name = "expression_value")

head(vst_final_D1)
nrow(vst_final_D1)

# Enhanced PCA plot
library(ggplot2)
library(pheatmap)
pca_plot_D1 <- plotPCA(vsdata_D1, intgroup = "Treatment") +
  theme_minimal() +
  labs(title = "PCA - Osimertinib 6w vs Untreated",
       subtitle = paste("n =", ncol(counts_D1), "samples")) +
  theme(plot.title = element_text(hjust = 0.5, size = 14),
        plot.subtitle = element_text(hjust = 0.5, size = 12))
pca_plot_D1

#saving PCA plot in pdf & png 
ggsave(("D1_PCA.pdf"), pca_plot_D1, width = 6, height = 5)
ggsave(("D1_PCA.png"), pca_plot_D1, width = 6, height = 5, dpi = 300)

# Sample distance heatmap
sample_distances_D1 <- dist(t(assay(vsdata_D1)))
sampleDistMatrix_D1 <- as.matrix(sample_distances_D1)
pheatmap::pheatmap(sampleDistMatrix_D1,
         clustering_distance_rows = sample_distances_D1,
         clustering_distance_cols = sample_distances_D1,
         main = "Sample-to-Sample Distances",
         filename = "SampleDistances_D1.png")
# Dispersion plot
plotDispEsts(dds_D1)

# Save as PDF
pdf("DispersionPlot_D1.pdf", width = 6, height = 5)
plotDispEsts(dds_D1)
dev.off()

# Or save as high-resolution PNG
png("DispersionPlot_D1.png", width = 2000, height = 1600, res = 300)
plotDispEsts(dds_D1)
dev.off()

# Get results with contrast
resultsNames(dds_D1)
res_D1 <- results(dds_D1, contrast = c("Treatment", "Osimertinib 6w", "Untreated"))
summary(res_D1)

# Apply LFC shrinkage
library(apeglm)
res_shrunk_D1 <- lfcShrink(dds_D1, coef=2, type="apeglm")

# Add gene annotations (ENSEMBL to SYMBOL)
library(dplyr)
library(AnnotationDbi)
library(org.Hs.eg.db)
gene_symbols_D1 <- mapIds(org.Hs.eg.db,
                      keys = rownames(res_shrunk_D1),
                      column = "SYMBOL",
                      keytype = "ENSEMBL",
                      multiVals = "first")

res_D1_annotated <- res_D1
res_D1_annotated$gene_symbol <- gene_symbols_D1  


as.data.frame(res_D1_annotated)

res_shrunk_annotated_D1 <- res_shrunk_D1
res_shrunk_annotated_D1$gene_symbol_D1 <- gene_symbols_D1

# Define significant genes
sig_overall_D1 <- res_shrunk_D1[which(res_shrunk_D1$padj < 0.05 & 
                                     abs(res_shrunk_D1$log2FoldChange) > 1), ]
head(sig_overall_D1)
# Separate up and down regulated genes
sig_genes_up_D1 <- rownames(res_shrunk_D1)[res_shrunk_D1$log2FoldChange > 1 & 
                                        res_shrunk_D1$padj < 0.05 & 
                                        !is.na(res_shrunk_D1$padj)]
# Extract gene symbols for your significant upregulated genes
gene_symbols_sig_up_D1 <- res_D1_annotated[sig_genes_up_D1, "gene_symbol"]

# Create a data frame with both gene IDs and symbols
sig_genes_up_D1_annotated <- data.frame(
  gene_id = sig_genes_up_D1,
  gene_symbol = gene_symbols_sig_up_D1
)
head(sig_genes_up_D1)
sig_genes_down_D1 <- rownames(res_shrunk_D1)[res_shrunk_D1$log2FoldChange < -1 & 
                                         res_shrunk_D1$padj < 0.05 & 
                                         !is.na(res_shrunk_D1$padj)]

print(paste("Significant DEGs:", nrow(sig_overall_D1)))
print(paste("Upregulated:", length(sig_genes_up_D1)))
print(paste("Downregulated:", length(sig_genes_down_D1)))

# MA Plot
plotMA(res_shrunk_D1, ylim=c(-5,5), main="Osimertinib vs Untreated")

pdf(("D1_MAplot.pdf"), width = 6, height = 5)
plotMA(res_shrunk_D1, ylim=c(-5,5), main="Osimertinib vs Untreated")
dev.off()

png(("D1_MAplot.png"), width = 2000, height = 1600, res = 300)
plotMA(res_shrunk_D1, ylim=c(-5,5), main="Osimertinib vs Untreated")
dev.off()

library(tibble)  # for rownames_to_column

# Select top 20 most significant genes for labeling
top_genes_data <- res_shrunk_annotated_D1 %>%
  as.data.frame() %>%
  rownames_to_column("gene_id") %>%
  filter(!is.na(padj) & !is.na(gene_symbol_D1)) %>%
  arrange(padj) %>%
  head(20)

top_gene_symbols_D1 <- top_genes_data$gene_symbol_D1

# Clean volcano plot
library(EnhancedVolcano)

volcano_plot_D1 <- EnhancedVolcano(res_shrunk_annotated_D1,
                lab = res_shrunk_annotated_D1$gene_symbol_D1,
                selectLab = top_gene_symbols_D1,
                x = 'log2FoldChange',
                y = 'padj',
                title = 'Osimertinib 6w vs Untreated',
                subtitle = paste('DEGs:', length(sig_genes_up_D1), 'up,', length(sig_genes_down_D1), 'down'),
                pCutoff = 0.05,
                FCcutoff = 1,
                pointSize = 1.5,
                labSize = 3.5,
                drawConnectors = TRUE,
                xlim = c(-6, 6))

volcano_plot_D1

ggsave(("D1_Volcano.pdf"), volcano_plot_D1, width = 6, height = 5)
ggsave(("D1_Volcano.png"), volcano_plot_D1, width = 6, height = 5, dpi = 300)

# Heatmap of top 50 significant genes
library(pheatmap)
top50_genes_D1 <- head(order(res_shrunk_annotated_D1$padj), 50)
png("Top50_Heatmap_D1.png", width = 2000, height = 1600, res = 300)
pheatmap::pheatmap(assay(vsdata_D1)[top50_genes_D1,], 
         annotation_col = as.data.frame(colData(dds_D1)[,"Treatment", drop=FALSE]), 
         scale = "row", 
         show_rownames = FALSE,
         main = "Top 50 Significant Genes",
         cluster_cols = TRUE)
dev.off()

library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)
library(dplyr)

# First, prepare top genes for heatmap
top_genes_heatmap_D1 <- res_shrunk_annotated_D1[!is.na(res_shrunk_annotated_D1$padj) & 
                                          res_shrunk_annotated_D1$padj < 0.05 & 
                                          res_shrunk_annotated_D1$baseMean > 50 & 
                                          abs(res_shrunk_annotated_D1$log2FoldChange) > 2,]

top_genes_heatmap_D1 <- head(as.data.frame(top_genes_heatmap_D1), 100)  # More genes for 115 samples
top_genes_heatmap_D1 <- top_genes_heatmap_D1[order(top_genes_heatmap_D1$log2FoldChange, decreasing = TRUE),]

# Get normalized count data from vst
mat_D1 <- assay(vsdata_D1)[rownames(top_genes_heatmap_D1), ]
mat_scaled_D1 <- t(apply(mat_D1, 1, scale))
colnames(mat_scaled_D1) <- colnames(mat_D1)

# Select top and bottom genes (increase for 115 samples)
num_keep <- 40  # Show more genes with 115 samples
rows_keep_D1 <- c(seq(1:num_keep), seq((nrow(mat_scaled_D1) - num_keep + 1), nrow(mat_scaled_D1)))

# Prepare annotation matrices
l2_val_D1 <- as.matrix(top_genes_heatmap_D1[rows_keep_D1,]$log2FoldChange)
colnames(l2_val_D1) <- "logFC"

mean_D1 <- as.matrix(top_genes_heatmap_D1[rows_keep_D1,]$baseMean)
colnames(mean_D1) <- "AveExpr"

# Set up color schemes
col_logFC_D1 <- colorRamp2(c(min(l2_val_D1), 0, max(l2_val_D1)), c("blue", "white", "red"))
col_AveExpr_D1 <- colorRamp2(c(quantile(mean_D1)[1], quantile(mean_D1)[4]), c("white", "red"))
col_zscore_D1 <- colorRamp2(c(-2, 0, 2), c("#4575b4", "white", "#d73027"))

# Add sample annotation for 115 samples
sample_annotation <- HeatmapAnnotation(
  Treatment = metadata_D1$Treatment,  # Assuming your metadata has treatment column
  col = list(Treatment = c("Untreated" = "gray", "Osimertinib 6w" = "darkgreen"))
)

# Main expression heatmap - adjusted for 115 samples
h1_D1 <- Heatmap(
  mat_scaled_D1[rows_keep_D1,], 
  cluster_rows = TRUE,
  cluster_columns = TRUE,               # keep clustering to see subgroups
  column_split = metadata_D1$Treatment, # enforce 2 main blocks (Untreated / Osimertinib 6w)
  top_annotation = sample_annotation,
  show_row_names = FALSE,
  show_column_names = FALSE,
  column_names_gp = gpar(fontsize = 6),
  column_title = "115 Samples",
  heatmap_legend_param = list(title_gp = gpar(fontsize = 10)),
  col = col_zscore_D1,
  name = "Z-score"
)

# LogFC heatmap
h2_D1 <- Heatmap(l2_val_D1, 
                 row_labels = top_genes_heatmap_D1$gene_symbol[rows_keep_D1],
                 cluster_rows = FALSE, 
                 name = "logFC", 
                 col = col_logFC_D1,
                 width = unit(15, "mm"),
                 row_names_gp = gpar(fontsize = 6),  # Smaller font for many genes
                 cell_fun = function(j, i, x, y, w, h, col) { 
                   grid.text(round(l2_val_D1[i, j], 1), x, y, gp = gpar(fontsize = 5))
                 })

# Average expression heatmap
h3_D1 <- Heatmap(mean_D1, 
                 row_labels = top_genes_heatmap_D1$gene_symbol[rows_keep_D1],
                 cluster_rows = FALSE, 
                 name = "AveExpr", 
                 col = col_AveExpr_D1,
                 width = unit(15, "mm"),
                 row_names_gp = gpar(fontsize = 6),
                 cell_fun = function(j, i, x, y, w, h, col) { 
                   grid.text(round(mean_D1[i, j], 0), x, y, gp = gpar(fontsize = 5))
                 })

# Combine heatmaps
h_D1 <- h1_D1 + h2_D1 + h3_D1
h_D1

pdf("AnnotatedHeatmap_D1.pdf", width = 10, height = 6)
draw(h_D1)
dev.off()

png("AnnotatedHeatmap_D1.png", width = 3000, height = 2000, res = 300)
draw(h_D1)
dev.off()

# Convert ENSEMBL to ENTREZ for KEGG analysis
library(DOSE)
library(enrichplot)

library(clusterProfiler)
entrez_genes_up_D1 <- mapIds(org.Hs.eg.db, sig_genes_up_D1, "ENTREZID", "ENSEMBL")
entrez_genes_down_D1 <- mapIds(org.Hs.eg.db, sig_genes_down_D1, "ENTREZID", "ENSEMBL")
entrez_genes_up_D1 <- entrez_genes_up_D1[!is.na(entrez_genes_up_D1)]
entrez_genes_down_D1 <- entrez_genes_down_D1[!is.na(entrez_genes_down_D1)]

# GO Analysis

go_up_D1 <- enrichGO(gene = sig_genes_up_D1,
                    OrgDb = org.Hs.eg.db,
                    keyType = "ENSEMBL",
                    ont = "ALL",  # BP, CC, MF
                    pAdjustMethod = "BH",
                    pvalueCutoff = 0.05)
top_go_up_D1 <- go_up_D1@result[1:10, ]
as.data.frame(top_go_up_D1)

go_down_D1 <- enrichGO(gene = sig_genes_down_D1,
                      OrgDb = org.Hs.eg.db,
                      keyType = "ENSEMBL",
                      ont = "ALL",
                      pAdjustMethod = "BH",
                      pvalueCutoff = 0.05)
top_go_down_D1 <- go_down_D1@result[1:10, ]
as.data.frame(top_go_down_D1)

# KEGG Pathway Analysis
kegg_up_D1 <- enrichKEGG(gene = entrez_genes_up_D1,
                         organism = 'hsa',
                         pvalueCutoff = 0.05,
                         pAdjustMethod = "BH")
as.data.frame(kegg_up_D1)
top_kegg_up_D1 <- kegg_up_D1@result[1:10, ]
as.data.frame(top_kegg_up_D1)

kegg_down_D1 <- enrichKEGG(gene = entrez_genes_down_D1,
                           organism = 'hsa',
                           pvalueCutoff = 0.05,
                           pAdjustMethod = "BH")
as.data.frame(kegg_down_D1)
top_kegg_down_D1 <- kegg_down_D1@result[1:10,]
as.data.frame(top_kegg_down_D1)

# Visualize pathway results
if(nrow(go_up_D1@result) > 0) {
  p1 <- barplot(go_up_D1, split="ONTOLOGY", showCategory = 8) + 
    facet_grid(ONTOLOGY~., scale="free") +
    labs(title = "Upregulated GO Terms")
  print(p1)
  ggsave("GO_Upregulated_D1.pdf", p1, width = 7, height = 6)
  ggsave("GO_Upregulated_D1.png", p1, width = 7, height = 6, dpi = 300)
}

if(nrow(go_down_D1@result) > 0) {
  p2 <- barplot(go_down_D1, split="ONTOLOGY", showCategory = 8) + 
    facet_grid(ONTOLOGY~., scale="free") +
    labs(title = "Downregulated GO Terms")
  print(p1)
  ggsave("GO_Downregulated_D1.pdf", p2, width = 7, height = 6)
  ggsave("GO_Downregulated_D1.png", p2, width = 7, height = 6, dpi = 300)
}


if(nrow(kegg_up_D1@result) > 0) {
  p3 <- dotplot(kegg_up_D1, showCategory = 15, title = "Upregulated KEGG Pathways")
  print(p2)
  ggsave("KEGG_Upregulated_D1.pdf", p3, width = 7, height = 6)
  ggsave("KEGG_Upregulated_D1.png", p3, width = 7, height = 6, dpi = 300)
}

if(nrow(kegg_down_D1@result) > 0) {
  p4 <- dotplot(kegg_down_D1, showCategory = 15, title = "Downregulated KEGG Pathways")
  print(p3)
  ggsave("KEGG_Downregulated_D1.pdf", p4, width = 7, height = 6)
  ggsave("KEGG_Downregulated_D1.png", p4, width = 7, height = 6, dpi = 300)
}

# Prepare ranked gene list
library(fgsea)
library(msigdbr)
gene_ranks_D1 <- res_shrunk_D1$log2FoldChange
names(gene_ranks_D1) <- rownames(res_shrunk_D1)
gene_ranks_D1 <- gene_ranks_D1[!is.na(gene_ranks_D1)]
gene_ranks_D1 <- sort(gene_ranks_D1, decreasing = TRUE)

# Get Hallmark pathways
hallmark_sets_D1 <- msigdbr(species = "Homo sapiens", category = "H")
hallmark_list_D1 <- split(hallmark_sets_D1$ensembl_gene, hallmark_sets_D1$gs_name)

# Run GSEA
fgsea_results_D1 <- fgsea(pathways = hallmark_list_D1,
                         stats = gene_ranks_D1,
                         minSize = 15,
                         maxSize = 500,
                         nperm = 1000)

#plotting top 3 pathways
print("Top 3 GSEA pathways:")
top_pathways_D1 <- head(fgsea_results_D1[order(pval), ], 3)
print(top_pathways_D1[, c("pathway", "pval", "NES")])

for(i in 1:3) {
  pathway_names_D1 <- top_pathways_D1$pathway[i]
  p <- plotEnrichment(hallmark_list_D1[[pathway_names_D1]], gene_ranks_D1) + 
    labs(title = paste("Pathway", i, ":", pathway_names_D1))
  print(p)
  ggsave("Top3_GSEA_D1.pdf", p, width = 7, height = 6)
  ggsave("Top3_GSEA_D1.png", p, width = 7, height = 6, dpi = 300)
}

head(fgsea_results_D1[order(pval), ], 10)

# Summary statistics table
summary_stats_D1 <- data.frame(
  Metric = c("Total Samples", "Osimertinib Treated", "Untreated", 
            "Genes Analyzed", "Significant DEGs", "Upregulated", "Downregulated"),
  Value = c(ncol(counts_D1), 
           sum(metadata_D1$Treatment == "Osimertinib 6w"),
           sum(metadata_D1$Treatment == "Untreated"),
           nrow(res_shrunk_D1),
           nrow(sig_overall_D1),
           length(sig_genes_up_D1),
           length(sig_genes_down_D1))
)

print("Analysis Summary:")
print(summary_stats_D1)

# Top genes table
top_genes_table_D1 <- data.frame(
  Gene_Symbol = gene_symbols_D1[rownames(head(sig_overall_D1[order(sig_overall_D1$padj),], 20))],
  ENSEMBL_ID = rownames(head(sig_overall_D1[order(sig_overall_D1$padj),], 20)),
  Log2FC = round(head(sig_overall_D1[order(sig_overall_D1$padj),], 20)$log2FoldChange, 3),
  Padj = format(head(sig_overall_D1[order(sig_overall_D1$padj),], 20)$padj, scientific = TRUE, digits = 3),
  Direction = ifelse(head(sig_overall_D1[order(sig_overall_D1$padj),], 20)$log2FoldChange > 0, "Up", "Down")
)

print("Top 20 Significant Genes:")
print(top_genes_table_D1)

# Export main results
write.csv(as.data.frame(res_shrunk_annotated_D1), "D1_annotated_results.csv", row.names = TRUE)
write.csv(sig_overall_D1, "D1_significant_genes.csv", row.names = TRUE)
write.csv(summary_stats_D1, "D1_summary_statistics.csv", row.names = FALSE)
write.csv(top_genes_table_D1, "D1_top_genes.csv", row.names = FALSE)

# Export normalized data
write.csv(res_D1_annotated, "D1_results.csv")
write.csv(counts(dds_D1, normalized=TRUE), "D1_normalized_counts.csv", row.names = TRUE)
write.csv(vst_final_D1, "D1_vst_data.csv", row.names = TRUE)

# Export pathway results
if(nrow(go_up_D1@result) > 0) {
  write.csv(go_up_D1@result, "D1_GO_upregulated.csv", row.names = FALSE)
}
if(nrow(go_down_D1@result) > 0) {
  write.csv(go_down_D1@result, "D1_GO_downregulated.csv", row.names = FALSE)
}
if(nrow(kegg_up_D1@result) > 0) {
  write.csv(kegg_up_D1@result, "D1_KEGG_upregulated.csv", row.names = FALSE)
}
if(nrow(kegg_down_D1@result) > 0) {
  write.csv(kegg_down_D1@result, "D1_KEGG_downregulated.csv", row.names = FALSE)
}

# Remove the leadingEdge column (which contains lists)
fgsea_results_clean <- fgsea_results_D1 %>%
  select(-leadingEdge)

write.csv(fgsea_results_clean, "D1_GSEA_hallmarks.csv", row.names = FALSE)