#Loading your work directory and RNASeq counts and filtering the counts
setwd("//files.wustl.edu/Shares/DOM/ONC/Hirbe_Lab/Diana/UBR5 KO RNASeq/analysis/HOM VS WT_2")
Counts <- read.csv("Counts.csv")
# Remove duplicate rows from Counts
Counts <- Counts[!duplicated(Counts[, 1]), ]
rownames(Counts) <- Counts[, 1]
Counts<- Counts [, -1]
# Calculate row means
row_means <- rowMeans(Counts)
# Order genes by row means in descending order
ordered_counts <- Counts[rev(order(row_means)), ]
#alternative code to the above is (ordered_counts <- Counts[order(row_means, decreasing = TRUE), ])
# Filter out rows with row means less than 10
filtered_counts <- ordered_counts[rowMeans(ordered_counts) >= 10, ]
#save filtered data frame
write.csv(filtered_counts, "filtered_counts.csv")
#prepare metadata telling R the conditions (columns)
metadata <- data.frame(
sample_id = colnames(filtered_counts), # Assuming you have loaded the filtered expression data
condition = c(rep("UBR5 WT", 3), rep("UBR5 HOM", 3)), # Treatment conditions
replicate = c(1, 2, 3, 1, 2, 3) # Sample replicates
)
metadata$condition <- factor(metadata$condition, levels = c("UBR5 WT", "UBR5 HOM"))
#Load DESEQ2 for normalization
library(DESeq2)
#Use the DESeqDataSetFromMatrix function from DESeq2 to create a DESeqDataSet object
dds <- DESeqDataSetFromMatrix(countData = filtered_counts,
colData = metadata,
design = ~ condition)
#Normailzation of RNASeq data
#Perform normalization and estimation of dispersions: Use the DESeq() function to perform normalization and estimation of dispersions.
dds <- DESeq(dds)
results <- results(dds, alpha = 0.05)
DEGs <- subset(results, abs(log2FoldChange) > 1 & padj < 0.05)
#save the de_genes data frame
write.csv(DEGs, file = "DEG_HOM_VS_WT.csv",)
write.csv(results, file = "DeseqResults_HOM_VS_WT.csv",)
#Visualizing your differentially expressed genes
#create volcano plot
library(ggplot2)
# Add column to classify genes as DEG or not
results_df <- as.data.frame(results)
results_df$gene <- rownames(results_df)
results_df$threshold <- "Unchanged"
results_df$threshold[results_df$padj < 0.05 & abs(results_df$log2FoldChange) > 1] <- "DEG"
library(ggrepel)
# Volcano plot
ggplot(results_df, aes(x = log2FoldChange, y = -log10(padj), color = threshold)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("Unchanged" = "grey", "DEG" = "red")) +
theme_minimal(base_size = 14) +
labs(title = "Volcano Plot: HOM vs WT",
x = "log2 Fold Change (HOM vs WT)",
y = "-log10 Adjusted p-value",
color = "Gene status") +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black")
# Select top 50 significant genes by padj
top100 <- results_df[order(results_df$padj), ][1:100, ]
# Volcano plot
ggplot(results_df, aes(x = log2FoldChange, y = -log10(padj), color = threshold)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("Unchanged" = "grey", "DEG" = "red")) +
theme_minimal(base_size = 14) +
labs(title = "Volcano Plot: HOM vs WT",
x = "log2 Fold Change (HOM vs WT)",
y = "-log10 Adjusted p-value",
color = "Gene status") +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
geom_text_repel(data = top100,
aes(label = gene),
size = 3,
max.overlaps = Inf,
box.padding = 0.3,
point.padding = 0.2,
segment.color = "grey50")
write.csv(top100, file = "Top100_HOM_vs_WT.csv", row.names = FALSE)
#Querrying genes of interest i.e UBR5 substrates as determined form the papers you have read
# Define genes of interest
genes_of_interest <- c("Egfr", "Hsp90ab1", "Map2k2", "Cerk", "Pdgfra", "Tyk2", "Jak1",
"Yap1", "Taz", "Kdr", "Aurka", "Pten", "Csf1r","Ptch1", "Smo", "Gli2", "Gli3", "Wnt10a", "Rac2", "Rspo2", "Apc",
"Cd274", "Pdcd1", "Id1", "Id3", "Cdh1", "Cdc73", "Hrpt2","Csf1","Golph3", "Cdk1", "Acsl4", "Ptk2b", "Akt1", "Akt2", "Akt3", "Pik3ca", "Pik3c2a", "Pik3cb" , "Pik3c3", "Pik3c2b", "Pik3cd", "Atmin", "Cdkn1a", "Cdk9", "Rela", "Nfkb1", "Nfkb2", "Capza1", "Stat1", "Stat3", "Irf1", "Irf3")
# Subset DEGs for these genes (case-sensitive match!)
genes_subset <- results_df[rownames(results_df) %in% genes_of_interest, ]
# Save genes of interest with stats
write.csv(genes_subset, file = "GenesOfInterest_HOM_vs_WT.csv", row.names = TRUE)
ggplot(results_df, aes(x = log2FoldChange, y = -log10(padj), color = threshold)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("Unchanged" = "grey", "DEG" = "red")) +
theme_minimal(base_size = 14) +
labs(title = "Volcano Plot: HET vs WT",
x = "log2 Fold Change (HET vs WT)",
y = "-log10 Adjusted p-value",
color = "Gene status") +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
geom_text_repel(
data = genes_subset,
aes(label = gene),
color = "black", # <-- force label text to black
size = 3,
max.overlaps = Inf,
box.padding = 0.3,
point.padding = 0.2,
segment.color = "grey50"
)
# Subset DEGs only
DEGs_df <- as.data.frame(DEGs)
DEGs_df$gene <- rownames(DEGs_df)
# Find overlap between DEGs and genes of interest
genes_subset <- DEGs_df[rownames(DEGs_df) %in% genes_of_interest, ]
# Save overlapping genes with stats
write.csv(genes_subset, file = "GenesOfInterest_DEGs_HOM_vs_WT.csv", row.names = TRUE)
# Volcano plot with labels ONLY for genes of interest that are DEGs
ggplot(results_df, aes(x = log2FoldChange, y = -log10(padj), color = threshold)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("Unchanged" = "grey", "DEG" = "red")) +
theme_minimal(base_size = 14) +
labs(title = "Volcano Plot: HET vs WT",
x = "log2 Fold Change (HET vs WT)",
y = "-log10 Adjusted p-value",
color = "Gene status") +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
geom_text_repel(
data = genes_subset,
aes(label = gene),
color = "black", # labels in black
size = 5,
max.overlaps = Inf,
box.padding = 0.3,
point.padding = 0.2,
segment.color = "grey50"
)
#Running pathway enrichment analysis to determine pathways enriched following UBR5 KO
if (!requireNamespace("clusterProfiler", quietly = TRUE)) {
BiocManager::install("clusterProfiler")
}
if (!requireNamespace("msigdbr", quietly = TRUE)) {
install.packages("msigdbr")
}
library(clusterProfiler)
library(msigdbr)
# Convert results to dataframe
res_df <- as.data.frame(results)
# Remove NA log2FC
res_df <- res_df[!is.na(res_df$log2FoldChange), ]
# Create named vector: names = gene symbols, values = log2FC
gene_list <- res_df$log2FoldChange
names(gene_list) <- rownames(res_df)
# Sort decreasing for clusterProfiler
gene_list <- sort(gene_list, decreasing = TRUE)
# Mouse Hallmark gene sets
hallmark_sets <- msigdbr(species = "Mus musculus", category = "H") # H = Hallmark
# Use as two-column dataframe: gs_name (pathway), gene_symbol
term2gene <- hallmark_sets[, c("gs_name", "gene_symbol")]
# Make sure your DESeq2 results have no NA log2FC
res_df <- as.data.frame(results)
res_df <- res_df[!is.na(res_df$log2FoldChange), ]
# Named vector: names = gene symbols, values = log2FC
gene_list <- res_df$log2FoldChange
names(gene_list) <- rownames(res_df)
gene_list <- sort(gene_list, decreasing = TRUE)
gsea_res <- GSEA(
geneList = gene_list,
TERM2GENE = term2gene, # <- must be dataframe, not list
pvalueCutoff = 0.1,
verbose = FALSE
)
# View top pathways
head(as.data.frame(gsea_res))
# Save results
write.csv(as.data.frame(gsea_res), "GSEA_Hallmark_Mouse_HOM_vs_WT.csv", row.names = FALSE)
library(enrichplot)
# Convert GSEA results to dataframe
gsea_df <- as.data.frame(gsea_res)
# Suppose the top (and only) enriched pathway:
top_pathway <- gsea_df$ID[5] # or use $Description if you prefer
# Classic GSEA plot for the top pathway
gseaplot2(
gsea_res,
geneSetID = top_pathway, # pathway ID
title = gsea_df$Description[5], # nice descriptive title
color = "red"
)
#PLOT HALLMARK PATHWAYS
library(ggplot2)
# Convert GSEA results to dataframe
gsea_df <- as.data.frame(gsea_res)
# Order pathways by NES (normalized enrichment score)
gsea_df <- gsea_df[order(gsea_df$NES, decreasing = TRUE), ]
# Plot ALL enriched pathways
ggplot(gsea_df, aes(x = reorder(Description, NES), y = NES, fill = -log10(p.adjust))) +
geom_col() +
coord_flip() +
labs(
title = "GSEA: All Enriched Hallmark Pathways",
x = "Pathway",
y = "Normalized Enrichment Score (NES)",
fill = "-log10 adj p-value"
) +
theme_minimal(base_size = 14)
#plot GSEA KEGG, GO, REACTOME --------------------------
# Prepare ranked gene list with Entrez IDs for GSEA
# --------------------------
library(clusterProfiler)
library(org.Mm.eg.db)
library(ReactomePA)
# Convert gene symbols to Entrez IDs
entrez_map <- bitr(names(gene_list), fromType="SYMBOL", toType="ENTREZID", OrgDb=org.Mm.eg.db)
gene_list_df <- merge(entrez_map, data.frame(log2FC = gene_list), by.x="SYMBOL", by.y="row.names")
gene_list_df <- gene_list_df[!duplicated(gene_list_df$ENTREZID), ]
gene_list_named <- gene_list_df$log2FC
names(gene_list_named) <- gene_list_df$ENTREZID
# Sort decreasing for GSEA
gene_list_named <- sort(gene_list_named, decreasing = TRUE)
# --------------------------
# 1) GSEA: KEGG Pathways
# --------------------------
gsea_kegg <- gseKEGG(
geneList = gene_list_named,
organism = "mmu",
minGSSize = 10,
pvalueCutoff = 0.1,
verbose = TRUE
)
# Save KEGG GSEA results
write.csv(as.data.frame(gsea_kegg), "GSEA_KEGG_HOM_vs_WT.csv", row.names = FALSE)
# Top 30 KEGG pathways barplot
library(enrichplot)
library(ggplot2)
# Convert gseaResult to dataframe to see top pathways
gsea_df <- as.data.frame(gsea_kegg)
# Select top 30 pathways by NES or pvalue
top30 <- gsea_df[order(gsea_df$NES, decreasing = TRUE)[1:30], ]
# Ridgeplot (shows enrichment distribution for multiple pathways)
ridgeplot(gsea_kegg, showCategory = 30) +
ggtitle("GSEA: KEGG Top 30 Pathways") +
theme_minimal(base_size = 14)
# Optional: classic GSEA plot for the top pathway
top_pathway <- top30$ID[1]
gseaplot2(gsea_kegg, geneSetID = top_pathway,
title = top30$Description[1], color = "red")
# --------------------------
# 2) GSEA: GO Biological Process (BP)
# --------------------------
gsea_go_bp <- gseGO(
geneList = gene_list_named,
OrgDb = org.Mm.eg.db,
ont = "ALL",
keyType = "ENTREZID",
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 0.1,
verbose = TRUE
)
# Save GO BP GSEA results
write.csv(as.data.frame(gsea_go_bp), "GSEA_GO_BP_HOM_vs_WT.csv", row.names = FALSE)
# Top 30 GO BP pathways barplot
barplot(gsea_go_bp, showCategory = 30, title = "GSEA: GO BP Top 30 Pathways")
# --------------------------
# 3) GSEA: Reactome Pathways
# --------------------------
gsea_reactome <- gsePathway(
geneList = gene_list_named,
organism = "mouse",
minGSSize = 10,
pvalueCutoff = 0.1,
verbose = TRUE
)
# Save Reactome GSEA results
write.csv(as.data.frame(gsea_reactome), "GSEA_Reactome_HOM_vs_WT.csv", row.names = FALSE)
# Top 30 Reactome pathways barplot
barplot(gsea_reactome, showCategory = 30, title = "GSEA: Reactome Top 30 Pathways")
# --------------------------
library(ggplot2)
library(dplyr)
# Convert GSEA Reactome results to dataframe
gsea_reactome_df <- as.data.frame(gsea_reactome)
# Select top 30 pathways by NES magnitude
top30_reactome <- gsea_reactome_df %>%
arrange(desc(abs(NES))) %>%
slice(1:30)
# Reorder for plotting (highest NES on top)
top30_reactome$Description <- factor(top30_reactome$Description, levels = rev(top30_reactome$Description))
# Plot barplot: NES on x-axis, pathways on y-axis, fill by -log10(padj)
ggplot(top30_reactome, aes(x = NES, y = Description, fill = -log10(p.adjust))) +
geom_bar(stat = "identity") +
scale_fill_gradient(low = "red", high = "darkred") +
theme_minimal(base_size = 14) +
labs(title = "GSEA: Top 30 Reactome Pathways",
x = "Normalized Enrichment Score (NES)",
y = "",
fill = "-log10(adj.p)")
#Querrying the xepression of genes within a particular pathway you may be interested in i.e type I interferon signaling
ifna_genes <- c(
"Ifi208","Mndal","Pdcd5-ps","Gstp-ps","Ifi203-ps","Rrp8","Ifi206","Chchd10","Inca1","Ppif",
"Topors","Rrn3","Eaf2","Ticam1","Unc5b","Bmyc","Cth","Pttg1ip","Selenos","Fgb","Raf1",
"Ggct","Tmbim6","Kcnq3","Scn2a","D1Pas1","Acvr1","Pawr","Acvr1b","Adora2a","Parp1","Parp2",
"Agt","Agtr2","Akt1","Aldh2","Alox12","Ivns1abp","Slc25a4","Slc25a5","Anxa6","Apaf1","App",
"Ar","Atf2","Atf3","Atf4","Atm","Atp2a1","Atp7a","Atp5if1","Avp","Bad","Bak1","Bax","Bcl10",
"Bcl2","Bcl2a1a","Bcl2a1b","Bcl2a1c","Bcl2a1d","Bcl2l1","Bcl2l10","Bcl2l2","Bcl3","Bdkrb2",
"Bdnf","Bex3","Bid","Hrk","Bik","Bcl2l11","Bmi1","Bmp4","Bmpr1b","Bnip3","Bnip3l","Brca1",
"Brca2","Birc6","Bub1","Hyou1","Camk2b","Casp12","Casp2","Casp3","Casp6","Casp8","Casp9",
"Ctnna1","Ctnnb1","Cav1","Runx3","Rb1cc1","Cck","Cd24a","Cd28","Cd3e","Cd44","Cd5","Cdk11b",
"Cdkn1a","Cdkn2d","Cebpb","Cflar","Clu","Ackr3","Col2a1","Cradd","Creb3","Crh","Crip1",
"Csf2","Csnk2a1","Ctsc","Ctsh","Cttn","Cx3cr1","Cycs","Cyct","Cyp1b1","Dab2","Dapk2",
"Dapk3","Daxx","Dbh","Ddit3","Ddx3x","Ddx5","E2f1","Ei24","Eif2ak3","Eno1","Epha2","Epo",
"Erbb3","Ercc2","Ptprv","Esr2","Eya1","Eya2","Eya3","Eya4","Fadd","Faf1","Fas","Fasl",
"Fcgr2b","Fem1b","Fga","Fgf10","Fgf2","Fgfr1","Fgfr2","Fgfr3","Fhit","Fxn","Tlr3","Fyn",
"Fzd1","Fzd9","G0s2","Gas1","Gata1","Gata4","Usp15","Gcg","Gdnf","Gclc","Gclm","Gnai2","Gnai3",
"Rack1","Gpx1","Pdia3","Gstp2","Gstp1","Hdac2","Htt","Hells","Hgf","Hic1","Hif1a","Hint1",
"Hipk1","Hipk2","Hmox1","Hnrnpk","Hras","Dnaja1","Hspb1","Hspa1b","Hyal2","Icam1","Ier3",
"Ifi203","Ifi204","Ifnb1","Ifng","Igf1","Cd74","Ikbkg","Il10","Il12a","Il18","Il1a","Il1b",
"Il2","Il3","Il4","Il7","Inhba","Inhbb","Ins2","Itga6","Itgav","Itpr1","Jak2","Jak3","Jun",
"Kcnq2","Klf4","Krt18","Krt8","Lck","Lcn2","Lgals3","Lmna","Lta","Ltb","Ltbr","Mfn2","Sgk3",
"Bbc3","Rtkn2","Smad3","Smad4","Mal","Bmf","Maz","Mbd4","Mcl1","Mdm2","Mdm4","Melk","Kitl",
"Mif","Mknk1","Mknk2","Mlh1","Mmp2","Mmp9","Mnt","Meis3","Msh2","Msh6","Msx1","Mapt","Muc1",
"Myc","Nck1","Nck2","Nf1","Nfe2l2","Ngf","Ngfr","Nkx3-1","Nodal","Nog","Nrp1","Nr4a2","Osm",
"Mybbp1a","P2rx4","P2rx7","P4hb","Igbp1","Pdk2","Pdpk1","Pdx1","Pea15a","Pik3r1","Prkca",
"Prkcd","Serpine1","Plaur","Pml","Pmp22","Pnp","Septin4","Polb","Pou4f1","Pou4f2","Ppard",
"Ppef2","Ppp1ca","Prkdc","Mapk8ip1","Prodh","Psen1","Psen2","Psme3","Pten","Ptgis","Ptgs2",
"Ptpn1","Ptpn2","Ptprc","Rad9a","Nlrp1a","Rb1","Rela","Ret","Ripk1","Uri1","Rnf7","Rock2",
"Rpl26","Rps7","S100a8","S100a9","Scg2","Cx3cl1","Cxcl12","Sfrp2","Spi1","Sfrp1","Sgpl1",
"Shh","Siah1a","Siah1b","Siah2","Skil","Snai2","Siglec1","Snai1","Sod1","Sod2","Sort1","Sp100",
"Spn","Spop","Src","Stk11","Pycr1","Stx4a","Trp53bp2","Syk","Nrg1","Tifab","Taf6","Tcf7l2",
"Rhot2","Hip1","Agap2","Prdx2","Flcn","Arrb2","Tmc8","Tert","G2e3","Ifi27l2b","Tgfb1","Tgfb2",
"Tgfbr1","Thbs1","Tlr4","Tlr6","Ccar2","Tnf","Tnfaip3","Tnfrsf10b","Tnfrsf1a","Tnfrsf1b",
"Cd27","Tnfsf11","Tnfsf12","Dedd","Cd40lg","Cd70","Ppp1r13b","Tpd52l1","Traf1","Traf2",
"Tnfsf10","Plscr1","Trp53","Trp63","Trp73","Tpt1","Tnfrsf4","Tnfsf4","Ubb","Umod","Kdm6a",
"Stk24","Vdac2","Vdr","Vegfa","Dap","Vhl","Vnn1","Mrtfa","Senp1","Wfs1","Wnt1","Pak2",
"Wnt4","Wnt5a","Xbp1","Traf7","Bag6","Gstp3","Xpa","Yap1","Zfp13","Pcgf2","Ifi207","Stradb",
"Pdk1","Madd","Trib3","Eif2a","Tmem161a","Usp28","Ifi209","Nox1","Il20ra","Atad5","Dido1",
"Faim","Map2k5","Mapk7","Prkra","Peli3","Rbck1","Zfp385b","Pak5","E2f2","Nanos3","Eda2r",
"AY074887","Map2k1","Map2k4","Map3k5","Map3k7","Mapk8","Mapk9","Creb3l1","Ppia","Casp8ap2",
"Ern2","Aifm1","Acvr1c","Ppp2r5c","Ell3","Nherf1","Serinc3","Rps3","Bcap31","Adcy10","Tnfrsf12a",
"Phlda3","Nbn","Cep63","Bag3","Zfp385a","Hip1r","Siva1","Ifnz","Ercc6","Tmem117","Tnfsf15",
"Ep300","Il19","Fnip2","Card9","Tmem102","Parl","Rrm2b","Gfral","Itprip","Eno1b","Acsl5",
"Mettl21c","Hdac1","Gsdma3","Ero1a","Fbxw7","Fbh1","Prkn","Chek2","Tnfsf14","Pdcd7","Ppp2r1a",
"Srpx","Bok","Zfp622","Acaa2","Ifi27","Atp2a3","Ube2k","Pla2g6","Psmd10","Nono","Asah2","Ifi214",
"Pde3a","Sh3glb1","Plagl2","Gsdme","Sfn","Lgals12","Ubqln1","Becn1","Stk3","Higd1a","Nupr1",
"Aatf","Pdcd5","Pdcd10","Mtch2","Ybx3","Foxo3","Gabarap","Ikbke","Ripk3","Gsk3b","Ankrd2",
"Mllt11","Park7","Marchf7","Noc2l","Jmy","Pidd1","Stk4","Pmaip1","Pias4","Sh3rf1","Rhot1",
"Stk25","Fignl1","Mapk8ip2","Gsk3a","Ifi213","Faiml","Nlrp1b","Ube4b","Perp","Moap1","Herpud1",
"Itm2c","Htra2","Zfp110","Arl6ip5","Txndc12","Ghitm","Eef1e1","Grina","Ing5","Snw1","Fis1",
"Pam16","Ptpmt1","Prelid1","Zmynd11","Timm50","Diablo","Cdip1","Lrrk2","Gskip","Bcl2l14",
"Pycard","Rnf186","Dele1","Dnajc10","Shisa5","Ndufa13","Armc10","Rffl","Dedd2","Erp29",
"Rnf41","Ddx47","Rps27l","Nacc2","Trap1","Coa8","Aen","Ndufs3","Mul1","Steap3","Tmem109",
"Ppm1f","Pink1","Zfas1","Zdhhc3","Chac1","Triap1","Fcmr","Dyrk2","Qrich1","Ing2","Dab2ip",
"Dapk1","Tmbim1","Tfpt","Fbxo7","Trim32","Fam162a","Plscr3","Bag5","Sfpq","Tmem238l","Tradd",
"Zswim2","Faim2","Rps6kb1","Uaca","Bclaf1","Nfatc4","Slc25a31","Bloc1s2","Ppp2r1b","Bbln",
"Dnm1l","Ddias","Syvn1","Opa1","Cyld","Wdr35","Ddit4","Pik3cb","Slc35f6","Usp47","Nme5",
"Tmem14a","Mff","Bcl2l12","Brsk2","Rnf183","Knl1","Styxl1","Dapl1","Gper1","Ifi27l2a",
"Il33","Nol3","Ern1","Tnfrsf23","Tnfrsf22","Trim39","Wwox","Rnf34","Selenok","Clca3a2",
"Nfkbiz","Sgpp1","Trem2","Trps1","Phip","Mpv17l","Wnt16","Sirt1","Tm2d1","Maged1","Hmgb2",
"Qars1","Deptor","Mael","Fgg","Kdm1a"
)
# Get normalized counts from DESeq2
norm_counts <- counts(dds, normalized=TRUE)
# Subset for IFN-α response genes (keep only genes present in your dataset)
ifna_counts <- norm_counts[rownames(norm_counts) %in% ifna_genes, ]
# Optionally, z-score normalize each gene for heatmap visualization
ifna_counts_z <- t(scale(t(ifna_counts)))
# Create annotation for columns
ann_col <- data.frame(
Condition = metadata$condition
)
rownames(ann_col) <- metadata$sample_id
library(pheatmap)
# Make sure columns are in your desired order
desired_order <- c("sample.WT_1", "sample.WT_2", "sample.WT_3",
"sample.neg.neg_1", "sample.neg.neg_2", "sample.neg.neg_3")
ifna_counts_z <- ifna_counts_z[, desired_order]
# Make sure annotation matches
ann_col <- ann_col[desired_order, , drop = FALSE]
# Heatmap
pheatmap(ifna_counts_z,
annotation_col = ann_col,
show_rownames = TRUE,
show_colnames = TRUE,
cluster_rows = TRUE,
cluster_cols = FALSE, # keep the column order fixed
scale = "row",
fontsize_row = 8,
main = "HALLMARK_APOPTOSIS")
# Subset DEGs for IFNa-related genes
ifna_DEGs <- DEGs_df[DEGs_df$gene %in% ifna_genes, ]
# Save to CSV
write.csv(ifna_DEGs, "GOBP_APOPTOSIS_HOM_vs_WT DEG.csv", row.names = FALSE)
# Save subset expression to CSV
write.csv(as.data.frame(ifna_counts),
file = "GOBP_APOPTOSIS_HOM VS WT.csv",
row.names = TRUE)
# Quick check
print(ifna_DEGs)
library(ggplot2)
library(ggrepel)
ggplot(results_df, aes(x = log2FoldChange, y = -log10(padj), color = threshold)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("Unchanged" = "grey", "DEG" = "red")) +
theme_minimal(base_size = 14) +
labs(
title = "HALLMARK_APOPTOSIS: HOM vs WT",
x = "log2 Fold Change (HOM vs WT)",
y = "-log10 Adjusted p-value",
color = "Gene status"
) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
# Label only IFNa DEGs
geom_text_repel(
data = ifna_DEGs,
aes(x = log2FoldChange, y = -log10(padj), label = gene),
inherit.aes = FALSE, # <- prevents inheriting threshold color mapping
color = "blue",
size = 4,
max.overlaps = Inf,
box.padding = 0.3,
point.padding = 0.2,
segment.color = "black"
)
#Heatmap for DEGs GOBP Apoptosis
# Subset DEGs for IFNa-related genes
ifna_DEGs <- DEGs_df[DEGs_df$gene %in% ifna_genes, ]
# Subset normalized counts to only DEGs
ifna_counts_DEG <- ifna_counts[rownames(ifna_counts) %in% ifna_DEGs$gene, ]
# Z-score normalize each DEG for heatmap
ifna_counts_DEG_z <- t(scale(t(ifna_counts_DEG)))
# Make sure columns are in your desired order
ifna_counts_DEG_z <- ifna_counts_DEG_z[, desired_order]
# Annotation matches columns
ann_col_DEG <- ann_col[desired_order, , drop = FALSE]
# Heatmap of only IFNa DEGs
pheatmap(ifna_counts_DEG_z,
annotation_col = ann_col_DEG,
show_rownames = TRUE,
show_colnames = TRUE,
cluster_rows = TRUE,
cluster_cols = FALSE, # keep column order fixed
scale = "row",
fontsize_row = 8,
main = "GOBP Apoptosis_DEGs Heatmap")
