Dataset 4

library(DESeq2)
library(ggplot2)

counts_D4 <- read.delim("counts_D4", header = TRUE, sep = "\t", row.names = 1)

total_counts_per_gene_D4 <- rowSums(counts_D4)
hist(total_counts_per_gene_D4, 
     breaks=100, 
     main="Distribution of Total Counts per Gene", 
     xlab="Total counts across all samples",
     ylab="Number of genes")
summary(total_counts_per_gene_D4)

#ignoring low expressed gene to remove noise & get better results

counts_D4 <- counts_D4[which(rowSums(counts_D4) > 50), ]

#calling in the metadata
metadata_D4 <- read.csv("metadata_D4.csv", header = TRUE)

counts_D4 <- round(counts_D4)
dds_D4 <- DESeqDataSetFromMatrix(countData = counts_D4,
                                colData = metadata_D4,        
                                design = ~ PHENOTYPE)
dds_D4

dds_D4 <- DESeq(dds_D4)
dds_D4

res_D4 <- results(dds_D4)
res_D4

#moderate filtering 
keep_D4 <- rowSums(counts(dds_D4) >= 20) >= 3
dds_D4 <- dds_D4[keep_D4,]

cat("Genes retained with moderate filtering:", nrow(dds_D4), "\n")

#Checking data quality after DESeq2
dds_D4 <- DESeq(dds_D4, minReplicatesForReplace = Inf)

resultsNames(dds_D4)
summary_res_D4 <- results(dds_D4, name = "PHENOTYPE_resistant_vs_parental")
summary(summary_res_D4)
plotDispEsts(dds_D4, main = "Dispersion Estimates: EBC-1 vs EBC-CR")

pdf("DispersionPlot_D4.pdf", width = 6, height = 5)
plotDispEsts(dds_D4)
dev.off()

png("DispersionPlot_D4.png", width = 2000, height = 1600, res = 300)
plotDispEsts(dds_D4)
dev.off()



vst_D4 <- vst(dds_D4, blind = FALSE)

vst_matrix_D4 <- assay(vst_D4)

head(vst_matrix_D4[1:5, 1:3])

library(reshape2)

vst_df_D4 <- as.data.frame(vst_matrix_D4)
vst_df_D4$gene_id <- rownames(vst_df_D4)

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

head(vst_final_D4)

We increased the counts threshold because of the imbalance in groups, ie. 1 vs 3. We needed to clear the noise for better results.We got good results.

pca_plot_D4 <- plotPCA(vst_D4, intgroup = "PHENOTYPE") +
  ggtitle("PCA: EBC-1 Parental vs EBC-CR Resistant (1 vs 3 Imbalance)") +
  theme_minimal() +
  theme(legend.position = "bottom") 
pca_plot_D4

ggsave(("D4_PCA.pdf"), pca_plot_D4, width = 6, height = 5)

ggsave(("D4_PCA.png"), pca_plot_D4, width = 6, height = 5, dpi = 300)

The PCA “Spread” is NOT Due to Poor Quality:

• Correlations measure overall gene expression similarity (Pearson correlation).
• PCA captures the dimensions of greatest variation (even if small in absolute terms).
• High correlations (0.87-0.97) with clear PCA separation = high-quality data with real biological differences.

resultsNames(dds_D4)

[1] "Intercept"                       "PHENOTYPE_resistant_vs_parental"

results_D4 <- results(dds_D4,
                     alpha = 0.01,           
                     lfcThreshold = 1,      
                     cooksCutoff = FALSE,    
                     contrast = c("PHENOTYPE", "resistant", "parental"))

print("Results summary with strict thresholds:")

[1] "Results summary with strict thresholds:"

summary(results_D4)

out of 12601 with nonzero total read count
adjusted p-value < 0.01
LFC > 1.00 (up)    : 141, 1.1%
LFC < -1.00 (down) : 9, 0.071%
outliers [1]       : 0, 0%
low counts [2]     : 978, 7.8%
(mean count < 38)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

results_D4

log2 fold change (MLE): PHENOTYPE resistant vs parental 
Wald test p-value: PHENOTYPE resistant vs parental 
DataFrame with 12601 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat    pvalue      padj
                <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000142655   42.6408      -0.910796  0.922302 -0.987525  0.557669  0.999995
ENSG00000149527   75.4277      -1.811092  0.742950 -2.437705  0.137556  0.999995
ENSG00000171621   68.4871      -1.451062  0.678824 -2.137612  0.253346  0.999995
ENSG00000173614  210.2283       0.129475  0.668358  0.193720  0.949146  0.999995
ENSG00000171729  325.1116      -0.920575  0.524908 -1.753784  0.560262  0.999995
...                   ...            ...       ...       ...       ...       ...
ENSG00000273730   34.6519       0.769113  0.835253  0.920814  0.625975        NA
ENSG00000301761   85.0955       0.704629  0.655156  1.075513  0.678584  0.999995
ENSG00000304745   19.9733      -0.181435  0.940890 -0.192833  0.912468        NA
ENSG00000276345  180.2394      -0.888292  0.680519 -1.305316  0.567956  0.999995
ENSG00000271254  114.7438      -0.471812  0.608137 -0.775832  0.815205  0.999995

tryCatch({
  results_shrunk_D4 <- lfcShrink(dds_D4, 
                                coef="PHENOTYPE_resistant_vs_parental", 
                                type="apeglm")
  print("LFC shrinkage successful")
}, error = function(e) {
  print("LFC shrinkage failed due to sample imbalance - using unshrunken results")
  results_shrunk_D4 <- results_D4
})

[1] "LFC shrinkage successful"

library(AnnotationDbi)
library(org.Hs.eg.db)

# Converting ENSEMBL IDs to gene symbols
gene_symbols_D4 <- mapIds(org.Hs.eg.db,
                          keys = rownames(dds_D4),
                          column = "SYMBOL",
                          keytype = "ENSEMBL",
                          multiVals = "first")

res_annotated_D4 <- results_D4
res_annotated_D4$gene_symbol <- gene_symbols_D4
as.data.frame(res_annotated_D4)

# Adding annotations to results
res_annotated_D4 <- results_shrunk_D4
res_annotated_D4$gene_symbol <- gene_symbols_D4

# Using strict criteria for significance
sig_genes_D4 <- res_annotated_D4[which(res_annotated_D4$padj < 0.01 & 
                                       abs(res_annotated_D4$log2FoldChange) > 1), ]

nrow(sig_genes_D4)

[1] 377

# Get top genes
top_genes_D4 <- head(sig_genes_D4[order(sig_genes_D4$padj), ], 20)

# Since we have 3 resistant samples, check consistency
normalized_counts_D4 <- counts(dds_D4, normalized = TRUE)

# Calculate fold change for each resistant sample vs parental
parental_counts <- normalized_counts_D4[, metadata_D4$PHENOTYPE == "parental"]
resistant_samples <- colnames(normalized_counts_D4)[metadata_D4$PHENOTYPE == "resistant"]

# Check consistency across resistant samples
fc_matrix <- sapply(resistant_samples, function(sample) {
  log2((normalized_counts_D4[, sample] + 1) / (parental_counts + 1))
})

# Find genes consistently up or down in all resistant samples
consistent_up <- rownames(fc_matrix)[apply(fc_matrix > 0.5, 1, all)]
consistent_down <- rownames(fc_matrix)[apply(fc_matrix < -0.5, 1, all)]

# Robust genes: significant AND consistent
sig_genes_up_D4 <- rownames(sig_genes_D4)[sig_genes_D4$log2FoldChange > 0]
sig_genes_down_D4 <- rownames(sig_genes_D4)[sig_genes_D4$log2FoldChange < 0]

robust_up_D4 <- intersect(sig_genes_up_D4, consistent_up)
robust_down_D4 <- intersect(sig_genes_down_D4, consistent_down)

cat("Significant genes:", nrow(sig_genes_D4), "\n")

Significant genes: 377 

cat("Robust upregulated (significant + consistent):", length(robust_up_D4), "\n")

Robust upregulated (significant + consistent): 326 

cat("Robust downregulated (significant + consistent):", length(robust_down_D4), "\n")

Robust downregulated (significant + consistent): 51 

# MA Plot
plotMA(results_shrunk_D4, ylim = c(-5, 5), main = "Resistant vs Parental")

pdf(("MA_D4.pdf"), width = 6, height = 5)

plotMA(results_shrunk_D4, ylim=c(-5,5), main="Resistant vs Parental")
dev.off()

png 
  2 

png(("MA_D4.png"), width = 2000, height = 1600, res = 300)
plotMA(results_shrunk_D4, ylim=c(-5,5), main="Resistant vs Parental")

dev.off()

png 
  2 

# Volcano plot
library(EnhancedVolcano)
volcano_plot <- EnhancedVolcano(res_annotated_D4, lab = res_annotated_D4$gene_symbol, 
                x = 'log2FoldChange', y = 'padj',
                title = 'EBC-CR Resistant vs EBC-1 Parental',
                pCutoff = 0.01,
                FCcutoff = 1)
ggsave(("D4_Volcano.pdf"), volcano_plot, width = 6, height = 5)

ggsave(("D4_Volcano.png"), volcano_plot, width = 6, height = 5, dpi = 300)

# Heatmap of significant genes
library(pheatmap)
if(nrow(sig_genes_D4) > 0) {
  top_sig_genes <- head(rownames(sig_genes_D4[order(sig_genes_D4$padj), ]), 50)
  pheatmap::pheatmap(assay(vst_D4)[top_sig_genes,], 
           annotation_col = as.data.frame(colData(dds_D4)[, "PHENOTYPE", drop = FALSE]), 
           scale = "row", 
           show_rownames = FALSE,
           main = "Top 50 Differentially Expressed Genes",
           filename = "DEGs_D4.png")
}

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

# Preparing top genes for enhanced heatmap - focus on robust genes
if(length(c(robust_up_D4, robust_down_D4)) > 10) {
  # Combine robust up and down genes
  robust_genes_D4 <- c(robust_up_D4, robust_down_D4)
  robust_genes_data <- res_annotated_D4[robust_genes_D4,]
  
  # Sort by fold change
  robust_genes_data <- robust_genes_data[order(robust_genes_data$log2FoldChange, decreasing = TRUE),]
  top_robust_D4 <- head(rownames(robust_genes_data), 50)
  
  # Getting normalized count data
  mat_D4 <- assay(vst_D4)[top_robust_D4, ]
  mat_scaled_D4 <- t(apply(mat_D4, 1, scale))
  colnames(mat_scaled_D4) <- colnames(mat_D4)
  
  # Selecting top and bottom genes
  num_keep <- min(25, floor(nrow(mat_scaled_D4)/2))
  rows_keep_D4 <- c(seq(1:num_keep), seq((nrow(mat_scaled_D4) - num_keep + 1), nrow(mat_scaled_D4)))
  
  # Preparing annotation matrices
  l2_val_D4 <- as.matrix(robust_genes_data[rows_keep_D4,]$log2FoldChange)
  colnames(l2_val_D4) <- "logFC"
  
  mean_D4 <- as.matrix(robust_genes_data[rows_keep_D4,]$baseMean)
  colnames(mean_D4) <- "AveExpr"
  
  # Setting up color schemes
  col_logFC_D4 <- colorRamp2(c(min(l2_val_D4), 0, max(l2_val_D4)), c("blue", "white", "red"))
  col_AveExpr_D4 <- colorRamp2(c(quantile(mean_D4)[1], quantile(mean_D4)[4]), c("white", "red"))
  col_zscore_D4 <- colorRamp2(c(-2, 0, 2), c("#4575b4", "white", "#d73027"))
  
  # Creating heatmaps
  h1_D4 <- Heatmap(mat_scaled_D4[rows_keep_D4,], 
                   cluster_rows = TRUE,
                   column_labels = colnames(mat_scaled_D4), 
                   name = "Z-score",
                   col = col_zscore_D4,
                   cluster_columns = TRUE,
                   show_row_names = FALSE,
                   column_names_gp = gpar(fontsize = 9),
                   column_title = "Robust Resistance Genes")
  
  h2_D4 <- Heatmap(l2_val_D4, 
                   row_labels = robust_genes_data$gene_symbol[rows_keep_D4],
                   cluster_rows = FALSE, 
                   name = "logFC", 
                   col = col_logFC_D4,
                   width = unit(15, "mm"),
                   show_column_names = TRUE,
                   row_names_gp = gpar(fontsize = 6),
                   column_names_gp = gpar(fontsize = 9),
                   cell_fun = function(j, i, x, y, w, h, col) { 
                     grid.text(round(l2_val_D4[i, j], 1), x, y, gp = gpar(fontsize = 7))
                   })
  
  h3_D4 <- Heatmap(mean_D4, 
                   row_labels = robust_genes_data$gene_symbol[rows_keep_D4],
                   cluster_rows = FALSE, 
                   name = "AveExpr", 
                   col = col_AveExpr_D4,
                   width = unit(15, "mm"),
                   show_row_names = TRUE,
                   show_column_names = TRUE,
                   column_names_gp = gpar(fontsize = 9),
                   row_names_gp = gpar(fontsize = 6),
                   cell_fun = function(j, i, x, y, w, h, col) { 
                     grid.text(round(mean_D4[i, j], 0), x, y, gp = gpar(fontsize = 7))
                   })
  
  # Combining all heatmaps
  h_D4 <- h1_D4 + h2_D4 + h3_D4
  h_D4
  
pdf("AnnotatedHeatmap_D4.pdf", width = 10, height = 6)
draw(h_D4)
dev.off()

png("AnnotatedHeatmap_D4.png", width = 3000, height = 2000, res = 300)
draw(h_D4)
dev.off()
}

library(DOSE)
library(enrichplot)
library(clusterProfiler)

# Using robust genes for pathway analysis
# Converting ENSEMBL to ENTREZ for KEGG analysis
entrez_genes_up_D4 <- mapIds(org.Hs.eg.db, robust_up_D4, "ENTREZID", "ENSEMBL")
entrez_genes_down_D4 <- mapIds(org.Hs.eg.db, robust_down_D4, "ENTREZID", "ENSEMBL")
entrez_genes_up_D4 <- entrez_genes_up_D4[!is.na(entrez_genes_up_D4)]
entrez_genes_down_D4 <- entrez_genes_down_D4[!is.na(entrez_genes_down_D4)]

# Running GO Analysis
  go_up <- enrichGO(gene = robust_up_D4,
                       OrgDb = org.Hs.eg.db,
                       keyType = "ENSEMBL",
                       ont = "ALL",
                       pAdjustMethod = "BH",
                       pvalueCutoff = 0.05)
  as.data.frame(go_up)
  top_go_up_D4 <- go_up@result[1:10, ]
as.data.frame(top_go_up_D4)
  
  go_up_D4 <- plot(barplot(go_up, showCategory = 15) + labs(title = "Upregulated GO Terms (Resistant vs Parental)"))
  ggsave("GO_up_D4.pdf", go_up_D4, width = 7, height = 6)
  ggsave("GO_up_D4.png", go_up_D4, width = 7, height = 6, dpi = 300)

  go_down <- enrichGO(gene = robust_down_D4,
                         OrgDb = org.Hs.eg.db,
                         keyType = "ENSEMBL",
                         ont = "ALL",
                         pAdjustMethod = "BH",
                         pvalueCutoff = 0.05)
  as.data.frame(go_down)
  go_down_D4 <- plot(barplot(go_down, showCategory = 15) + labs(title = "Downregulated GO Terms (Resistant vs Parental)"))
  ggsave("GO_down_D4.pdf", go_down_D4, width = 7, height = 6)
  ggsave("GO_down_D4.png", go_down_D4, width = 7, height = 6, dpi = 300)

# Running KEGG Pathway Analysis
  kegg_up <- enrichKEGG(gene = entrez_genes_up_D4,
                           organism = 'hsa',
                           pvalueCutoff = 0.05,
                           pAdjustMethod = "BH")
  as.data.frame(kegg_up)
  top_kegg_up_D4 <- kegg_up@result[1:10, ]
as.data.frame(top_kegg_up_D4)
  
  
 kegg_up_D4 <-  plot(barplot(kegg_up, showCategory = 15) + labs(title = "Upregulated KEGG Pathways (Resistant vs Parental)"))
  ggsave("KEGG_up_D4.pdf", kegg_up_D4, width = 7, height = 6)
  ggsave("KEGG_up_D4.png", kegg_up_D4, width = 7, height = 6, dpi = 300)


# Add this after the KEGG analysis attempt
  kegg_down <- enrichKEGG(gene = entrez_genes_down_D4,
                             organism = 'hsa',
                             pvalueCutoff = 0.1,
                             pAdjustMethod = "BH")
  as.data.frame(kegg_down)
  
  if(nrow(as.data.frame(kegg_down)) > 0) {
    print(as.data.frame(kegg_down))
    plot(barplot(kegg_down, showCategory = 15) + 
         labs(title = "Downregulated KEGG Pathways (Resistant vs Parental)"))
  } else {
    cat("No KEGG pathways enriched for downregulated genes (n =", 
        length(entrez_genes_down_D4), "genes)\n")
  }

library(fgsea)
library(msigdbr)

# Preparing ranked gene list for GSEA
gene_ranks_D4 <- results_shrunk_D4$log2FoldChange
names(gene_ranks_D4) <- rownames(results_shrunk_D4)
gene_ranks_D4 <- gene_ranks_D4[!is.na(gene_ranks_D4)]
gene_ranks_D4 <- sort(gene_ranks_D4, decreasing = TRUE)

# Getting Hallmark pathways for GSEA
hallmark_sets_D4 <- msigdbr(species = "Homo sapiens", category = "H")
hallmark_list_D4 <- split(hallmark_sets_D4$ensembl_gene, hallmark_sets_D4$gs_name)

# Running GSEA analysis
fgsea_results_D4 <- fgsea(pathways = hallmark_list_D4,
                          stats = gene_ranks_D4,
                          minSize = 15,
                          maxSize = 500)

# Viewing top results
print("Top GSEA Results (Resistant vs Parental):")
print(head(fgsea_results_D4[order(pval), ], 10))

# Plotting top 3 most significant pathways
print("Top 3 GSEA pathways:")
top_pathways_D4 <- head(fgsea_results_D4[order(pval), ], 3)
print(top_pathways_D4[, c("pathway", "pval", "NES")])

for(i in 1:3) {
  pathway_name <- top_pathways_D4$pathway[i]
  pathways_D4 <- plotEnrichment(hallmark_list_D4[[pathway_name]], gene_ranks_D4) + 
    labs(title = paste("Pathway", i, ":", pathway_name))
  print(pathways_D4)
  ggsave("Top3_GSEA_D4.pdf", pathways_D4, width = 7, height = 6)
  ggsave("Top3_GSEA_D4.png", pathways_D4, width = 7, height = 6, dpi = 300)
}

# Summary statistics table for D4
summary_stats_D4 <- data.frame(
  Metric = c("Total Samples", "Genes Analyzed", "Significant DEGs", "Upregulated genes", "Downregulated genes"),
  Value = c(ncol(counts_D4), 
            nrow(res_D4),
            nrow(sig_genes_D4),
            length(sig_genes_up_D4),
            length(sig_genes_down_D4))
)

summary_stats_D4

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

print("Top 20 Significant Genes")

[1] "Top 20 Significant Genes"

print(top_genes_table_D4)

NA

# Export main results
write.csv(as.data.frame(res_annotated_D4), "D4_annotated_results.csv", row.names = TRUE)
write.csv(sig_genes_D4, "D4_sig_genes.csv", row.names = TRUE)
write.csv(summary_stats_D4, "D4_summary_statistics.csv", row.names = FALSE)
write.csv(top_genes_table_D4, "D4_top20_genes.csv", row.names = FALSE)

# Export normalized data
write.csv(counts(dds_D4, normalized=TRUE), "D4_normalized_counts.csv", row.names = TRUE)
write.csv(assay(vst_D4), "D4_vst_data.csv", row.names = TRUE)

# Export pathway results
  write.csv(go_up@result, "D4_GO_upregulated.csv", row.names = FALSE)

  write.csv(go_down@result, "D4_GO_downregulated.csv", row.names = FALSE)


  write.csv(kegg_up@result, "D4_KEGG_upregulated.csv", row.names = FALSE)

# Remove the leadingEdge column (which contains lists)
fgsea_results_clean_D4 <- fgsea_results_D4 %>%
  select(-leadingEdge)

write.csv(fgsea_results_clean_D4, "D4_GSEA_hallmarks.csv", row.names = FALSE)