Differential Gene Analysis and Visualization for HSC Methylome Data

A Gene Ontology (GO) enrichment analysis was conducted on differential genes derived from both male and female conditions. The datasets were imported and preprocessed, followed by the application of the clusterProfiler package for enrichment analysis. Subsequently, significant GO terms were identified and visualized in a tailored bubble plot, with colors signifying different ontologies. The finalized results were documented in both EXCEL and PNG formats.

install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("enrichplot")
if(!requireNamespace("org.Mm.eg.db", quietly = TRUE)) {
  BiocManager::install("org.Mm.eg.db")
}

# Load necessary libraries
library(clusterProfiler)
library(enrichplot)
library(org.Mm.eg.db)
library(AnnotationDbi)

# Read differential genes and GO terms
male_diff_genes <- read.csv("male_n_reworked.csv", header = TRUE, stringsAsFactors = FALSE)
#go_terms_list <- read.csv("male_n_reworked.csv", header = TRUE, stringsAsFactors = FALSE)
# Assuming the genes are in a column named 'genes' and GO terms are in a column named 'go_terms'
male_gene_vector_clean <- na.omit(male_diff_genes$GeneSymbol)

# # Convert to EntrezIDs
# male_entrez_ids <- select(org.Mm.eg.db, keys = male_gene_vector_clean, keytype = "SYMBOL", columns = "ENTREZID")
# # Filter out any NAs or genes not found
# male_entrez_ids_clean <- na.omit(male_entrez_ids$ENTREZID)

# Perform GO enrichment analysis
# Using 'org.Hs.eg.db' as an example for human. Replace with appropriate organism db if different.
male_ego <- enrichGO(gene = male_gene_vector_clean, 
                    OrgDb = org.Mm.eg.db, 
                    keyType = "SYMBOL",
                    ont = "ALL", 
                    pAdjustMethod = "none",
                    pvalueCutoff = 1.0,
                    qvalueCutoff = 1.0, 
                    readable = TRUE)                    
print(male_ego$ID)
View(male_ego)
write.csv(male_ego, "male_n_reworked_GO.csv", row.names = FALSE)

#NOTE that we should manually delete the second line of csv-file since the line is empty in the original Excel-file.
# Read differential genes and GO terms
female_diff_genes <- read.csv("female_n_reworked.csv", header = TRUE, stringsAsFactors = FALSE)
#go_terms_list <- read.csv("female_n_reworked.csv", header = TRUE, stringsAsFactors = FALSE)
# Assuming the genes are in a column named 'genes' and GO terms are in a column named 'go_terms'
female_gene_vector_clean <- na.omit(female_diff_genes$GeneSymbol)

# Perform GO enrichment analysis
# Using 'org.Hs.eg.db' as an example for human. Replace with appropriate organism db if different.
female_ego <- enrichGO(gene = female_gene_vector_clean, 
                    OrgDb = org.Mm.eg.db, 
                    keyType = "SYMBOL",
                    ont = "ALL", 
                    pAdjustMethod = "none",
                    pvalueCutoff = 1.0,
                    qvalueCutoff = 1.0, 
                    readable = TRUE)
print(female_ego$ID)
View(female_ego)
write.csv(female_ego, "female_n_reworked_GO.csv", row.names = FALSE)

#~/Tools/csv2xls-0.4/csv_to_xls.py male_n_reworked_GO.csv female_n_reworked_GO.csv -d',' -o GO_enrichments.xls


# ---- calculate the significant gene-list ----
male_ego_ <- enrichGO(gene = male_gene_vector_clean, 
                    OrgDb = org.Mm.eg.db, 
                    keyType = "SYMBOL",
                    ont = "ALL", 
                    pAdjustMethod = "BH",
                    readable = TRUE)
female_ego_ <- enrichGO(gene = female_gene_vector_clean, 
                    OrgDb = org.Mm.eg.db, 
                    keyType = "SYMBOL",
                    ont = "ALL", 
                    pAdjustMethod = "BH",
                    readable = TRUE)
# Merge and deduplicate
merged_description <- unique(c(male_ego_$Description, female_ego_$Description))
merged_ID <- unique(c(male_ego_$ID, female_ego_$ID))

#go_terms_list <- read.csv("go_terms.csv", header = TRUE, stringsAsFactors = FALSE)
#go_terms_vector <- go_terms_list$go_terms
#go_terms_vector <- c("GO:0090140", "GO:0006664", "GO:0046488", "GO:1903509")
# Filter the results by the provided GO terms
male_filtered_ego <- subset(male_ego, male_ego$ID %in% merged_ID)
female_filtered_ego <- subset(female_ego, female_ego$ID %in% merged_ID)
female_filtered_ego$Sex <- "Female"
male_filtered_ego$Sex <- "Male"
mydat <- rbind(male_filtered_ego, female_filtered_ego)

library(ggplot2)
library(dplyr)
library(magrittr)
library(tidyr)
library(forcats)

#mydat <- read.csv2("GO_BP_enrichment_genebody.txt", sep="\t", header=TRUE)
#mydat$FoldEnrichment <- as.numeric(mydat$FoldEnrichment)
#mydat$GeneRatio <- sapply(as.character(mydat$GeneRatio), function(x) eval(parse(text=x)))
#mydat$p.adjust <- as.numeric(as.character(mydat$FDR))
mydat$Sex <- factor(mydat$Sex, levels=c("Male","Female")) 
mydat$Description <- factor(mydat$Description, levels=rev(merged_description))
mydat$ONTOLOGY <- factor(mydat$ONTOLOGY, levels=c("BP","CC", "MF"))
mydat$log_p.adjust <- -log10(mydat$p.adjust)

# Define a custom scale function and its inverse
rescale_fun <- function(x) {
  return(-log10(x))
}

inverse_fun <- function(x) {
  return(10^(-x))
}

# Use the custom scale function in the ggplot
png("bubble.png", width=800, height=860)
ggplot(mydat, aes(y = Description, x = Sex, size = p.adjust)) + 
geom_point(aes(color = ONTOLOGY), alpha = 1.0) + 
labs(x = "", y = "") + 
theme(axis.text.x = element_text(angle = 30, vjust = 0.5)) + 
theme(axis.text = element_text(size = 20)) + 
theme(legend.text = element_text(size = 20)) + 
theme(legend.title = element_text(size = 20)) + 
scale_size_continuous(trans = scales::trans_new(name = "custom", transform = rescale_fun, inverse = inverse_fun),
                      breaks = c(1, 0.5, 0.1, 0.05, 0.01, 0.001), 
                      labels = c("1", "0.5", "0.1", "0.05", "0.01", "0.001"),
                      range = c(1, 13)) + 
guides(color = guide_legend(override.aes = list(size = 10))) +
scale_color_manual(values = c("BP" = "blue", "CC" = "red", "MF" = "green"))
dev.off()

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默