---
title: "GO_cluster"
author: "Luise A. Seeker"
date: "08/11/2021"
output: html_document
---
# introduction
I tested GSEA and gene ontology both for molecular functions and biological
pathways and different visualisations. In general GO works better for my data
than GSEA so therefore in later versions of this document, there may be no
code for GSEA any more.
```{r}
library(here)
library(tibble)
library(Seurat)
library(dplyr)
library(ggplot2)
library(tidyr)
#GO & GSEA
library(biomaRt)
library(org.Hs.eg.db)
library(ReactomePA)
library(msigdbr)
library(fgsea)
library(NMF)
library(ggupset)
library(clusterProfiler)
library(enrichplot)
```
# Read in data
```{r}
nad_ol <- readRDS(here("data",
"single_nuc_data",
"oligodendroglia",
"srt_oligos_and_opcs_LS.RDS"))
Idents(nad_ol) <- "ol_clusters_named"
opcs <- subset(nad_ol, idents = c("OPC_A", "OPC_B"))
Idents(opcs) <- "Tissue"
csc_ba4_mark <- FindMarkers(opcs, ident.1 = "CSC",
ident.2 = "BA4",
only.pos = FALSE,
min.pct = 0.25,
test.use = "MAST")
```
# Prepare hallmarks for GSEA and gene sets
```{r, eval = FALSE}
hallmark <- msigdbr(species = "Homo sapiens", category = "H")
geneset <- hallmark %>% split(x = .$gene_symbol, f = .$gs_name)
```
```{r}
ba4 <- subset(csc_ba4_mark, csc_ba4_mark$avg_log2FC < 0)
csc <- subset(csc_ba4_mark, csc_ba4_mark$avg_log2FC >0)
# Prepare marts to convert to entrez
listMarts()
ensembl = useMart("ensembl", dataset="hsapiens_gene_ensembl",
host = "www.ensembl.org")
listDatasets(ensembl)
attributes = listAttributes(ensembl)
genes_ba4 <- rownames(ba4)
# convert to ENTREZID
genes_ent <- bitr(genes_ba4, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
# Tissue
ck_kegg_tissue <- compareCluster(geneCluster = genes_ent, fun = enrichKEGG)
ck_kegg_tissue <- setReadable(ck_kegg_tissue, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_tissue)
dotplot(ck_kegg_tissue)
```
```{r}
genes_csc <- rownames(csc)
# convert to ENTREZID
genes_ent <- bitr(genes_csc, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
# Tissue
ck_kegg_tissue <- compareCluster(geneCluster = genes_ent, fun = enrichKEGG)
ck_kegg_tissue <- setReadable(ck_kegg_tissue, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_tissue)
dotplot(ck_kegg_tissue)
```
```{r}
fil_data <- subset(csc, csc$p_val_adj < 0.05)
genes <- rownames(fil_data)
fil_data$abs_avg_log2FC <- abs(fil_data$avg_log2FC)
fil_data$gene <- rownames(fil_data)
# Prepare input data
ranks <- fil_data %>% dplyr::select(gene, abs_avg_log2FC)
ranks <- deframe(ranks)
#RUN GO using clusterProfiler
fil_data$Gene <- fil_data$gene
#remove any duplicates (sanity check for me)
genemodulesGO <- fil_data[!duplicated(fil_data$Gene),]
Biomart_gencode_ensembl84_biotypes <- getBM(attributes=c("ensembl_gene_id",
"entrezgene_id",
"gene_biotype",
"hgnc_symbol"),
filters = "hgnc_symbol",
values = fil_data$Gene,
mart = ensembl,
uniqueRows = FALSE)
Biomart_gencode_ensembl84_biotypes[, 'gene_biotype'] <-
as.factor(Biomart_gencode_ensembl84_biotypes[,'gene_biotype'])
#Filter genes
#biotype_all_dataset <- subset(Biomart_gencode_ensembl84_biotypes, hgnc_symbol
# %in% nad_ol@assays$RNA@var.features)
#entrezID <- subset(biotype_all_dataset, biotype_all_dataset$hgnc_symbol
# %in% nad_ol@assays$RNA@var.features)
# if(nrow(biotype_all_dataset) == 0){
biotype_all_dataset <- subset(Biomart_gencode_ensembl84_biotypes, hgnc_symbol
%in% rownames(opcs@assays$RNA))
entrezID <- subset(biotype_all_dataset, biotype_all_dataset$hgnc_symbol
%in% rownames(opcs@assays$RNA))
# }
entrezmatched <- entrezID[match(genemodulesGO$Gene,entrezID$hgnc_symbol),]
#you might need to remove NAs
entrezmatched <- entrezmatched[! apply(entrezmatched[,c(1,3)], 1,function (x) anyNA(x)),]
allLLIDs <- entrezmatched$entrezgene
modulesReactome <- enrichPathway(gene=as.character(allLLIDs),organism="human",
pvalueCutoff=0.01,
qvalueCutoff = 0.3,
pAdjustMethod = "none",
readable=T)
if(nrow(as.data.frame(modulesReactome)) > 0){
dot_plot <- dotplot(modulesReactome, showCategory=20)
pdf(file = here(go_dir, plot_label),
width = 12,
height = 10)
print(dot_plot)
dev.off()
x2 <- pairwise_termsim(modulesReactome)
emap_plot <- emapplot(x2)
pdf(file = here(emap_dir, plot_label),
width = 12,
height = 10)
print(emap_plot)
dev.off()
#Alternative GO clasissification
ent_uni <- rownames(nad_ol@assays$RNA)
ent_uni_entrez <- bitr(ent_uni, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
ent_uni <- ent_uni_entrez$ENTREZID
# Molecular function
go_mf <- enrichGO(gene = as.character(allLLIDs),
universe = ent_uni,
OrgDb = org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 1,
qvalueCutoff = 1,
readable = TRUE)
go_mf_plot <- dotplot(go_mf, showCategory = 20)
pdf(file = here(go_dir_mol_f, plot_label),
width = 12,
height = 10)
print(go_mf_plot)
dev.off()
# GO biological process
ego <- enrichGO(gene = as.character(allLLIDs),
OrgDb = org.Hs.eg.db,
ont = "BP",
universe = ent_uni )
#head(ego)
if(nrow(ego) > 1){
dp2 <- dotplot(ego, showCategory=20)
pdf(file = here(go_dot, plot_label),
width = 12,
height = 10)
print(dp2)
dev.off()
## convert gene ID to Symbol
edox <- setReadable(ego, 'org.Hs.eg.db', 'ENTREZID')
cnet_pl <- cnetplot(edox, foldChange=ranks, circular = TRUE,
colorEdge = TRUE)
pdf(file = here(go_dir_cnetpl, plot_label),
width = 14,
height = 12)
print(cnet_pl)
dev.off()
edox2 <- enrichplot::pairwise_termsim(edox)
if(nrow(edox2) > 1){
tree_plot <- treeplot(edox2)
pdf(file = here(go_dir_tree, plot_label), # The directory you want to save the file in
width = 16, # The width of the plot in inches
height = 10) # The height of the plot in inches
print(tree_plot)
dev.off()
}
u_plot <- upsetplot(ego)
pdf(file = here(go_dir_upset_plot, plot_label),
width = 15,
height = 10)
print(u_plot)
dev.off()
if(nrow(edox) > 1){
hm <- heatplot(edox, foldChange=ranks, showCategory=5)
pdf(file = here(go_dir_hm, plot_label), # The directory you want to save the file in
width = 15, # The width of the plot in inches
height = 10) # The height of the plot in inches
print(hm)
dev.off()
}
}
edox_mf <- setReadable(go_mf, 'org.Hs.eg.db', 'ENTREZID')
edox2_mf <- pairwise_termsim(edox_mf)
tree_plot_mf <- treeplot(edox2_mf)
pdf(file = here(go_dir_mol_f_tree, plot_label),
width = 16,
height = 10)
print(tree_plot_mf)
dev.off()
}
}
}
```
# for non oligos
```{r}
files <- c( here("outs",
"astrocytes",
"tissue_marker_list_astro",
"astro_tissue_marker.csv"),
here("outs",
"astrocytes",
"age_marker_list_astro",
"astro_age_marker.csv"),
here("outs",
"astrocytes",
"sex_marker_list_astro",
"astro_sex_marker.csv"),
here("outs",
"astrocytes",
"age_sex_marker_list_astro",
"astro_age_sex_marker.csv"),
###
here("outs",
"microglia_macrophages",
"tissue_marker_list_mm",
"mm_tissue_marker.csv"),
here("outs",
"microglia_macrophages",
"age_marker_list_mm",
"mm_age_marker.csv"),
here("outs",
"microglia_macrophages",
"sex_marker_list_mm",
"mm_sex_marker.csv"),
here("outs",
"microglia_macrophages",
"age_sex_marker_list_mm",
"mm_age_sex_marker.csv"),
###
here("outs",
"vascular_cells",
"tissue_marker_list_vasc",
"vasc_tissue_mark.csv"),
here("outs",
"vascular_cells",
"age_marker_list_vasc",
"vasc_age_marker.csv"),
here("outs",
"vascular_cells",
"sex_marker_list_vasc",
"vasc_sex_marker.csv"),
here("outs",
"vascular_cells",
"age_sex_marker_list_vasc",
"vasc_age_sex_mark.csv"),
###
here("outs",
"neurons",
"tissue_marker_list_neuron",
"neuron_tissue_marker.csv"),
here("outs",
"neurons",
"age_marker_list_neuron",
"neuron_age_marker.csv"),
here("outs",
"neurons",
"sex_marker_list_neuron",
"neuron_sex_marker.csv"),
here("outs",
"neurons",
"age_sex_marker_list_neuron",
"age_sex_neurons_marker.csv")
)
single_nuc_data <- c(here("data",
"single_nuc_data",
"astrocytes",
"HCA_astrocytes.RDS"),
here("data",
"single_nuc_data",
"microglia",
"HCA_microglia.RDS"),
here("data",
"single_nuc_data",
"vascular_cells",
"HCA_vascular_cells.RDS"),
here("data",
"single_nuc_data",
"neurons",
"HCA_neurons.RDS")
)
```
```{r}
named_list <- list()
for(i in 1: length(files)){
data <- read.csv(files[i])
cluster_lev <- levels(as.factor(data$cluster))
for(k in 1:length(cluster_lev)){
curr_data <- subset(data, data$cluster == cluster_lev[k])
cluster <- curr_data[["cluster"]][1]
cell_t <- sapply(strsplit(files[i], "/"), `[`, 12)
cell_t <- sapply(strsplit(cell_t, "_marker"), `[`, 1)
label <- paste0(cell_t, "_", cluster)
plot_label <- paste0(label, ".pdf")
fil_data <- subset(curr_data, curr_data$p_val_adj < 0.05 &
curr_data$avg_log2FC > 0)
genes <- fil_data$gene
# convert to ENTREZID
genes_ent <- bitr(genes, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
named_list[[k]] <- genes_ent$ENTREZID
names(named_list)[k] <- label
}
if(i == 1){
return_list <- named_list
}else{
return_list <- append(return_list, named_list)
}
named_list <- list()
}
str(return_list)
tissue_list_astro <- return_list[1:3]
age_list_astro <- return_list[4:5]
sex_list_astro <- return_list[6:7]
age_sex_list_astro <- return_list[8:11]
tissue_list_mm <- return_list[12:14]
age_list_mm <- return_list[15:16]
sex_list_mm <- return_list[17:18]
age_sex_list_mm <- return_list[19:22]
tissue_list_vasc <- return_list[23:25]
age_list_vasc <- return_list[26:27]
sex_list_vasc <- return_list[28:29]
age_sex_list_vasc <- return_list[30:33]
tissue_list_neur <- return_list[34:36]
age_list_neur <- return_list[37:38]
sex_list_neur <- return_list[39:40]
age_sex_list_neur <- return_list[41:44]
# ASTROCYTES
# Tissue
ck_kegg_tissue <- compareCluster(geneCluster = tissue_list_astro, fun = enrichKEGG)
ck_kegg_tissue <- setReadable(ck_kegg_tissue, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_tissue)
dotplot(ck_kegg_tissue)
# Age
ck_kegg_age <- compareCluster(geneCluster = age_list_astro, fun = enrichKEGG)
ck_kegg_age <- setReadable(ck_kegg_age, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age)
dotplot(ck_kegg_age)
# Age*sex
ck_kegg_age_sex <- compareCluster(geneCluster = age_sex_list_astro, fun = enrichKEGG)
ck_kegg_age_sex <- setReadable(ck_kegg_age_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age_sex)
dotplot(ck_kegg_age_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# sex
ck_kegg_sex <- compareCluster(geneCluster = sex_list_astro, fun = enrichKEGG)
ck_kegg_sex <- setReadable(ck_kegg_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_sex)
dotplot(ck_kegg_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# MICROGLIA
# Tissue
ck_kegg_tissue <- compareCluster(geneCluster = tissue_list_mm, fun = enrichKEGG)
ck_kegg_tissue <- setReadable(ck_kegg_tissue, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_tissue)
dotplot(ck_kegg_tissue)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# Age
ck_kegg_age <- compareCluster(geneCluster = age_list_mm, fun = enrichKEGG)
ck_kegg_age <- setReadable(ck_kegg_age, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age)
dotplot(ck_kegg_age)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# Age*sex
ck_kegg_age_sex <- compareCluster(geneCluster = age_sex_list_mm, fun = enrichKEGG)
ck_kegg_age_sex <- setReadable(ck_kegg_age_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age_sex)
dotplot(ck_kegg_age_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# sex
ck_kegg_sex <- compareCluster(geneCluster = sex_list_mm, fun = enrichKEGG)
ck_kegg_sex <- setReadable(ck_kegg_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_sex)
dotplot(ck_kegg_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# VASCULAR
# Tissue
ck_kegg_tissue <- compareCluster(geneCluster = tissue_list_vasc, fun = enrichKEGG)
ck_kegg_tissue <- setReadable(ck_kegg_tissue, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_tissue)
dotplot(ck_kegg_tissue)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# Age
ck_kegg_age <- compareCluster(geneCluster = age_list_vasc, fun = enrichKEGG)
ck_kegg_age <- setReadable(ck_kegg_age, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age)
dotplot(ck_kegg_age)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# Age*sex
ck_kegg_age_sex <- compareCluster(geneCluster = age_sex_list_vasc, fun = enrichKEGG)
ck_kegg_age_sex <- setReadable(ck_kegg_age_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age_sex)
dotplot(ck_kegg_age_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# sex
ck_kegg_sex <- compareCluster(geneCluster = sex_list_vasc, fun = enrichKEGG)
ck_kegg_sex <- setReadable(ck_kegg_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_sex)
dotplot(ck_kegg_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# neurons
# Tissue
ck_kegg_tissue <- compareCluster(geneCluster = tissue_list_neur, fun = enrichKEGG)
ck_kegg_tissue <- setReadable(ck_kegg_tissue, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_tissue)
dotplot(ck_kegg_tissue)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# Age
ck_kegg_age <- compareCluster(geneCluster = age_list_neur, fun = enrichKEGG)
ck_kegg_age <- setReadable(ck_kegg_age, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age)
dotplot(ck_kegg_age)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# Age*sex
ck_kegg_age_sex <- compareCluster(geneCluster = age_sex_list_neur, fun = enrichKEGG)
ck_kegg_age_sex <- setReadable(ck_kegg_age_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_age_sex)
dotplot(ck_kegg_age_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
# sex
ck_kegg_sex <- compareCluster(geneCluster = sex_list_neur, fun = enrichKEGG)
ck_kegg_sex <- setReadable(ck_kegg_sex, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck_kegg_sex)
dotplot(ck_kegg_sex) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
```
```{r}
for(i in 1: length(files)){
data <- read.csv(files[i])
cluster_lev <- levels(as.factor(data$cluster))
for(k in 1:length(cluster_lev)){
curr_data <- subset(data, data$cluster == cluster_lev[k])
cluster <- curr_data[["cluster"]][1]
cluster <- gsub(" ", "_", cluster)
cell_t <- sapply(strsplit(files[i], "/"), `[`, 12)
cell_t <- sapply(strsplit(cell_t, "_mark"), `[`, 1)
label <- paste0(cell_t, "_", cluster)
plot_label <- paste0(label, ".pdf")
cell_type <- sapply(strsplit(cell_t, "_"), `[`, 1)
fil_data <- subset(curr_data, curr_data$p_val_adj < 0.05 &
curr_data$avg_log2FC > 0)
if(nrow(fil_data) > 0){
genes <- fil_data$gene
# Prepare input data
ranks <- fil_data %>% dplyr::select(gene, avg_log2FC)
ranks <- deframe(ranks)
#RUN GO using clusterProfiler
fil_data$Gene <- fil_data$gene
#remove any duplicates (sanity check for me)
genemodulesGO <- fil_data[!duplicated(fil_data$Gene),]
Biomart_gencode_ensembl84_biotypes <- getBM(attributes=c("ensembl_gene_id",
"entrezgene_id",
"gene_biotype",
"hgnc_symbol"),
filters = "hgnc_symbol",
values = fil_data$Gene,
mart = ensembl,
uniqueRows = FALSE)
Biomart_gencode_ensembl84_biotypes[, 'gene_biotype'] <-
as.factor(Biomart_gencode_ensembl84_biotypes[,'gene_biotype'])
if(cell_type == "astro"){
srt <- readRDS(single_nuc_data[1])
}else if(cell_type == "mm"){
srt <- readRDS(single_nuc_data[2])
}else if(cell_type == "vasc"){
srt <- readRDS(single_nuc_data[3])
}else if(cell_type == "neuron" | cell_type == "age"){
srt <- readRDS(single_nuc_data[4])
}else{
print("check srt.")
}
#Filter genes
#biotype_all_dataset <- subset(Biomart_gencode_ensembl84_biotypes, hgnc_symbol
# %in% srt@assays$RNA@var.features)
#entrezID <- subset(biotype_all_dataset, biotype_all_dataset$hgnc_symbol
# %in% srt@assays$RNA@var.features)
# if(nrow(biotype_all_dataset) == 0){
biotype_all_dataset <- subset(Biomart_gencode_ensembl84_biotypes, hgnc_symbol
%in% rownames(srt@assays$RNA))
entrezID <- subset(biotype_all_dataset, biotype_all_dataset$hgnc_symbol
%in% rownames(srt@assays$RNA))
# }
entrezmatched <- entrezID[match(genemodulesGO$Gene,entrezID$hgnc_symbol),]
#you might need to remove NAs
entrezmatched <- entrezmatched[! apply(entrezmatched[,c(1,3)], 1,function (x) anyNA(x)),]
allLLIDs <- entrezmatched$entrezgene
modulesReactome <- enrichPathway(gene=as.character(allLLIDs),organism="human",
pvalueCutoff=0.01,
qvalueCutoff = 0.3,
pAdjustMethod = "none",
readable=T)
if(nrow(as.data.frame(modulesReactome)) > 0){
dot_plot <- dotplot(modulesReactome, showCategory=20)
pdf(file = here(go_dir, plot_label),
width = 12,
height = 10)
print(dot_plot)
dev.off()
x2 <- pairwise_termsim(modulesReactome)
emap_plot <- emapplot(x2)
pdf(file = here(emap_dir, plot_label),
width = 12,
height = 10)
print(emap_plot)
dev.off()
#Alternative GO clasissification
ent_uni <- rownames(srt@assays$RNA)
ent_uni_entrez <- bitr(ent_uni, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
ent_uni <- ent_uni_entrez$ENTREZID
# Molecular function
go_mf <- enrichGO(gene = as.character(allLLIDs),
universe = ent_uni,
OrgDb = org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 1,
qvalueCutoff = 1,
readable = TRUE)
go_mf_plot <- dotplot(go_mf, showCategory = 20)
pdf(file = here(go_dir_mol_f, plot_label),
width = 12,
height = 10)
print(go_mf_plot)
dev.off()
# GO biological process
ego <- enrichGO(gene = as.character(allLLIDs),
OrgDb = org.Hs.eg.db,
ont = "BP",
universe = ent_uni )
#head(ego)
if(nrow(ego) > 1){
dp2 <- dotplot(ego, showCategory=20)
pdf(file = here(go_dot, plot_label),
width = 12,
height = 10)
print(dp2)
dev.off()
## convert gene ID to Symbol
edox <- setReadable(ego, 'org.Hs.eg.db', 'ENTREZID')
cnet_pl <- cnetplot(edox, foldChange=ranks, circular = TRUE,
colorEdge = TRUE)
pdf(file = here(go_dir_cnetpl, plot_label),
width = 14,
height = 12)
print(cnet_pl)
dev.off()
# edox2 <- enrichplot::pairwise_termsim(edox)
# if(nrow(edox2) > 1){
# tree_plot <- treeplot(edox2)
#pdf(file = here(go_dir_tree, plot_label), # The directory you want to save the file in
# width = 16, # The width of the plot in inches
# height = 10) # The height of the plot in inches
# print(tree_plot)
#dev.off()
u_plot <- upsetplot(ego)
pdf(file = here(go_dir_upset_plot, plot_label),
width = 15,
height = 10)
print(u_plot)
dev.off()
}
if(nrow(edox) > 1){
hm <- heatplot(edox, foldChange=ranks, showCategory=5)
pdf(file = here(go_dir_hm, plot_label), # The directory you want to save the file in
width = 15, # The width of the plot in inches
height = 10) # The height of the plot in inches
print(hm)
dev.off()
}
}
edox_mf <- setReadable(go_mf, 'org.Hs.eg.db', 'ENTREZID')
edox2_mf <- pairwise_termsim(edox_mf)
tree_plot_mf <- treeplot(edox2_mf)
pdf(file = here(go_dir_mol_f_tree, plot_label),
width = 16,
height = 10)
print(tree_plot_mf)
dev.off()
}
}
}
```
# Create and save files that contain GO output for oligos
```{r}
oligo_files
srt <- nad_ol
plot_df <- data.frame()
for(i in 1: length(oligo_files)){
data <- read.csv(oligo_files[i])
#srt <- readRDS(single_nuc_data[i])
cluster_lev <- levels(as.factor(data$cluster))
for(k in 1:length(cluster_lev)){
curr_data <- subset(data, data$cluster == cluster_lev[k])
cluster <- curr_data[["cluster"]][1]
#cell_t <- sapply(strsplit(oligo_files[i], "/"), `[`, 10)
cell_t <- "oligos_"
plot_label <- paste0(cell_t, "cluster", "GO_term_age_sex_region.csv")
fil_data <- subset(curr_data, curr_data$p_val_adj < 0.05 &
curr_data$avg_log2FC > 0)
#RUN GO using clusterProfiler
fil_data$Gene <- fil_data$gene
#remove any duplicates (sanity check for me)
genemodulesGO <- fil_data[!duplicated(fil_data$Gene),]
Biomart_gencode_ensembl84_biotypes <- getBM(attributes=c("ensembl_gene_id",
"entrezgene_id",
"gene_biotype",
"hgnc_symbol"),
filters = "hgnc_symbol",
values = fil_data$Gene,
mart = ensembl,
uniqueRows = FALSE)
Biomart_gencode_ensembl84_biotypes[, 'gene_biotype'] <-
as.factor(Biomart_gencode_ensembl84_biotypes[,'gene_biotype'])
biotype_all_dataset <- subset(Biomart_gencode_ensembl84_biotypes, hgnc_symbol
%in% rownames(srt@assays$RNA))
entrezID <- subset(biotype_all_dataset, biotype_all_dataset$hgnc_symbol
%in% rownames(srt@assays$RNA))
entrezmatched <- entrezID[match(genemodulesGO$Gene,entrezID$hgnc_symbol),]
#you might need to remove NAs
entrezmatched <- entrezmatched[! apply(entrezmatched[,c(1,3)], 1,function (x) anyNA(x)),]
allLLIDs <- entrezmatched$entrezgene
allLLIDs<-allLLIDs[!is.na(allLLIDs)]
ent_uni <- rownames(nad_ol@assays$RNA)
ent_uni_entrez <- bitr(ent_uni, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
ent_uni <- ent_uni_entrez$ENTREZID
ego <- enrichGO(gene = as.character(allLLIDs),
OrgDb = org.Hs.eg.db,
ont = "BP",
universe = ent_uni )
ego_df <- as.data.frame(ego)
if(nrow(ego_df)> 0){
ego_df$Cell_type <- cell_t
ego_df$cluster <- cluster_lev[k]
if(nrow(plot_df) < 1){
plot_df <- ego_df
}else{
plot_df <- rbind(plot_df, ego_df)
}
}
}
}
dir.create(here("outs",
"GO",
"tissue_age_sex_table_output",
"oligos"), recursive = TRUE)
write.csv(plot_df, here("outs",
"GO",
"tissue_age_sex_table_output",
"oligos",
paste0(plot_label)))
```
# Create and save files that contain GO output for non_oligos
```{r}
plot_df <- data.frame()
for(i in 1: length(files)){
data <- read.csv(files[i])
cell_t <- sapply(strsplit(files[i], "/"), `[`, 10)
if(cell_t == "astrocytes"){
srt <- readRDS(single_nuc_data[1])
}else if(cell_t == "microglia_macrophages"){
srt <- readRDS(single_nuc_data[2])
}else if(cell_t == "vascular_cells"){
srt <- readRDS(single_nuc_data[3])
}else if(cell_t == "neurons"){
srt <- readRDS(single_nuc_data[4])
}else{
print("check srt.")
}
cluster_lev <- levels(as.factor(data$cluster))
for(k in 1:length(cluster_lev)){
curr_data <- subset(data, data$cluster == cluster_lev[k])
cluster <- curr_data[["cluster"]][1]
plot_label <- paste0(cell_t, "_", cluster, "_GO_term_age_sex_region.csv")
fil_data <- subset(curr_data, curr_data$p_val_adj < 0.05 &
curr_data$avg_log2FC > 0)
if(nrow(fil_data) > 0){
#RUN GO using clusterProfiler
fil_data$Gene <- fil_data$gene
#remove any duplicates (sanity check for me)
genemodulesGO <- fil_data[!duplicated(fil_data$Gene),]
Biomart_gencode_ensembl84_biotypes <- getBM(attributes=c("ensembl_gene_id",
"entrezgene_id",
"gene_biotype",
"hgnc_symbol"),
filters = "hgnc_symbol",
values = fil_data$Gene,
mart = ensembl,
uniqueRows = FALSE)
Biomart_gencode_ensembl84_biotypes[, 'gene_biotype'] <-
as.factor(Biomart_gencode_ensembl84_biotypes[,'gene_biotype'])
biotype_all_dataset <- subset(Biomart_gencode_ensembl84_biotypes, hgnc_symbol
%in% rownames(srt@assays$RNA))
entrezID <- subset(biotype_all_dataset, biotype_all_dataset$hgnc_symbol
%in% rownames(srt@assays$RNA))
entrezmatched <- entrezID[match(genemodulesGO$Gene,entrezID$hgnc_symbol),]
#you might need to remove NAs
entrezmatched <- entrezmatched[! apply(entrezmatched[,c(1,3)], 1,function (x) anyNA(x)),]
allLLIDs <- entrezmatched$entrezgene
allLLIDs<-allLLIDs[!is.na(allLLIDs)]
ent_uni <- rownames(nad_ol@assays$RNA)
ent_uni_entrez <- bitr(ent_uni, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
ent_uni <- ent_uni_entrez$ENTREZID
ego <- enrichGO(gene = as.character(allLLIDs),
OrgDb = org.Hs.eg.db,
ont = "BP",
universe = ent_uni )
ego_df <- as.data.frame(ego)
if(nrow(ego_df)> 0){
ego_df$Cell_type <- cell_t
ego_df$cluster <- cluster_lev[k]
if(nrow(plot_df) < 1){
plot_df <- ego_df
}else{
plot_df <- rbind(plot_df, ego_df)
}
}
}
}
}
dir.create(here("outs",
"GO",
"tissue_age_sex_table_output",
"non_oligos"), recursive = TRUE)
write.csv(plot_df, here("outs",
"GO",
"tissue_age_sex_table_output",
"non_oligos",
"non_oligos_GO_term_age_sex_region.csv"))
```
# Session Info
```{r}
sessionInfo()
```