Skip to content

GitLab

Commit 31eabd95 authored by Monah Abou Alezz's avatar Monah Abou Alezz
Browse files

update scRNAseq scripts

parent 99a0892e
Hide whitespace changes
Inline Side-by-side
Showing
with 169 additions and 267 deletions
+169 -267
scRNAseq_mixed/seurat_analysis.R
View file @ 31eabd95
...@@ -11,78 +11,21 @@ suppressPackageStartupMessages(library(ggrepel)) ...@@ -11,78 +11,21 @@ suppressPackageStartupMessages(library(ggrepel))
suppressPackageStartupMessages(library(clusterProfiler)) suppressPackageStartupMessages(library(clusterProfiler))
suppressPackageStartupMessages(library(RColorBrewer)) suppressPackageStartupMessages(library(RColorBrewer))
## load seurat data ## load seurat data
h48_final <- readRDS("data/mixed_h48_final.rds") tp <- c("h48","d4")
d4_final <- readRDS("data/mixed_D4_final.rds") for (t in tp) {
Idents(h48_final) <- "RNA_snn_h.orig.ident_res.1.2" sc_obj <- readRDS(paste0("data/mixed_",t,"_final.rds"))
Idents(d4_final) <- "RNA_snn_h.orig.ident_res.1.2" Idents(sc_obj) <- "Cell_Type"
h48_final@meta.data[['orig.ident']] <- recode_factor( sc_obj$Cell_type <- factor(
h48_final@meta.data[["orig.ident"]], sc_obj$Cell_Type,
"Sample3" = "AAV9", levels = c("Astrocytes",
"Sample5" = "Spk", "Neurons",
"Sample7" = "UT") "Oligodendrocytes",
h48_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]] <- recode_factor( "Immature.Astrocytes",
h48_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]], "Immature.Neurons",
"0" = "Undifferentiated", "OPC",
"1" = "Astrocytes", "Undifferentiated"))
"2" = "Astrocytes", assign(paste0(t,"_final"), sc_obj)
"3" = "Astrocytes", }
"4" = "OPC",
"5" = "Oligodendrocytes",
"6" = "Astrocytes",
"7" = "Astrocytes",
"8" = "Oligodendrocytes",
"9" = "Immature.Neurons",
"10" = "Oligodendrocytes",
"11" = "Immature.Neurons",
"12" = "Immature.Astrocytes",
"13" = "Immature.Astrocytes",
"14" = "Astrocytes",
"15" = "Neurons",
"16" = "Oligodendrocytes",
"17" = "OPC")
h48_final$RNA_snn_h.orig.ident_res.1.2 <- factor(
h48_final$RNA_snn_h.orig.ident_res.1.2,
levels = c("Astrocytes",
"Neurons",
"Oligodendrocytes",
"Immature.Astrocytes",
"Immature.Neurons",
"OPC",
"Undifferentiated"))
d4_final@meta.data[['orig.ident']] <- recode_factor(
d4_final@meta.data[["orig.ident"]],
"Sample2" = "AAV2",
"Sample4" = "AAV9",
"Sample6" = "Spk",
"Sample8" = "UT")
d4_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]] <- recode_factor(
d4_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]],
"0" = "Undifferentiated",
"1" = "Oligodendrocytes",
"2" = "Neurons",
"3" = "Astrocytes",
"4" = "Immature.Neurons",
"5" = "OPC",
"6" = "Astrocytes",
"7" = "Oligodendrocytes",
"8" = "Oligodendrocytes",
"9" = "Immature.Neurons",
"10" = "Oligodendrocytes",
"11" = "Oligodendrocytes",
"12" = "Astrocytes",
"13" = "Immature.Astrocytes",
"14" = "OPC",
"15" = "Astrocytes",
"16" = "Astrocytes")
d4_final$RNA_snn_h.orig.ident_res.1.2 <- factor(
d4_final$RNA_snn_h.orig.ident_res.1.2,
levels = c("Astrocytes",
"Neurons",
"Oligodendrocytes",
"Immature.Astrocytes",
"Immature.Neurons",
"OPC",
"Undifferentiated"))
ggplotColours <- function(n = 6, h = c(0, 360) + 15){ ggplotColours <- function(n = 6, h = c(0, 360) + 15){
if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65) hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
...@@ -90,7 +33,7 @@ ggplotColours <- function(n = 6, h = c(0, 360) + 15){ ...@@ -90,7 +33,7 @@ ggplotColours <- function(n = 6, h = c(0, 360) + 15){
## umap ## umap
DimPlot(object = h48_final, DimPlot(object = h48_final,
reduction = "umap.harmony.orig.ident", reduction = "umap.harmony.orig.ident",
group.by = "RNA_snn_h.orig.ident_res.1.2", group.by = "Cell_Type",
label = T, label = T,
repel = T, repel = T,
pt.size = 0.2, pt.size = 0.2,
...@@ -105,7 +48,7 @@ DimPlot(object = h48_final, ...@@ -105,7 +48,7 @@ DimPlot(object = h48_final,
values = ggplotColours(n = 7)) values = ggplotColours(n = 7))
## donut plot ## donut plot
ggplot_object_clusters_labels <- melt(table(h48_final$RNA_snn_h.orig.ident_res.1.2)) ggplot_object_clusters_labels <- melt(table(h48_final$Cell_Type))
colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count") colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count)) totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
totalB$percentage <- totalB$percentage*100 totalB$percentage <- totalB$percentage*100
...@@ -159,7 +102,7 @@ calc_helper <- function(object,genes){ ...@@ -159,7 +102,7 @@ calc_helper <- function(object,genes){
aav9 <- subset(h48_final, orig.ident=="AAV9") aav9 <- subset(h48_final, orig.ident=="AAV9")
gfp_percentage_cell_type <- PrctCellExpringGene(aav9, gfp_percentage_cell_type <- PrctCellExpringGene(aav9,
genes ="GFP", genes ="GFP",
group.by = "RNA_snn_h.orig.ident_res.1.2") group.by = "Cell_Type")
gfp_percentage_cell_type$Percentage <- gfp_percentage_cell_type$Cell_proportion*100 gfp_percentage_cell_type$Percentage <- gfp_percentage_cell_type$Cell_proportion*100
gfp_percentage_cell_type <- gfp_percentage_cell_type[,c(4,3)] gfp_percentage_cell_type <- gfp_percentage_cell_type[,c(4,3)]
colnames(gfp_percentage_cell_type) <- c("Percentage", "Cell_Types") colnames(gfp_percentage_cell_type) <- c("Percentage", "Cell_Types")
...@@ -201,7 +144,7 @@ h48_aav9_ut <- subset(h48_final, orig.ident == "AAV9" | orig.ident == "UT") ...@@ -201,7 +144,7 @@ h48_aav9_ut <- subset(h48_final, orig.ident == "AAV9" | orig.ident == "UT")
d4_aav9_ut <- subset(d4_final, orig.ident == "AAV9" | orig.ident == "UT") d4_aav9_ut <- subset(d4_final, orig.ident == "AAV9" | orig.ident == "UT")
DimPlot(object = h48_aav9_ut, DimPlot(object = h48_aav9_ut,
reduction = "umap.harmony.orig.ident", reduction = "umap.harmony.orig.ident",
group.by = "RNA_snn_h.orig.ident_res.1.2", group.by = "Cell_Type",
split.by = "orig.ident", split.by = "orig.ident",
label = F, label = F,
repel = T, repel = T,
...@@ -218,7 +161,7 @@ DimPlot(object = h48_aav9_ut, ...@@ -218,7 +161,7 @@ DimPlot(object = h48_aav9_ut,
values = ggplotColours(n = 7)) values = ggplotColours(n = 7))
DimPlot(object = d4_aav9_ut, DimPlot(object = d4_aav9_ut,
reduction = "umap.harmony.orig.ident", reduction = "umap.harmony.orig.ident",
group.by = "RNA_snn_h.orig.ident_res.1.2", group.by = "Cell_Type",
split.by = "orig.ident", split.by = "orig.ident",
label = F, label = F,
repel = T, repel = T,
...@@ -234,12 +177,11 @@ DimPlot(object = d4_aav9_ut, ...@@ -234,12 +177,11 @@ DimPlot(object = d4_aav9_ut,
"Undifferentiated"), "Undifferentiated"),
values = ggplotColours(n = 7)) values = ggplotColours(n = 7))
samples <- c("AAV9", "UT") samples <- c("AAV9", "UT")
tp <- c("h48","d4")
for (t in tp) { for (t in tp) {
obj_t <- get(paste0(t,"_final")) obj_t <- get(paste0(t,"_final"))
for (i in samples) { for (i in samples) {
obj_i <- subset(obj_t, orig.ident == i) obj_i <- subset(obj_t, orig.ident == i)
ggplot_object_clusters_labels <- melt(table(obj_i$RNA_snn_h.orig.ident_res.1.2)) ggplot_object_clusters_labels <- melt(table(obj_i$Cell_Type))
colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count") colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count)) totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
totalB$percentage <- totalB$percentage*100 totalB$percentage <- totalB$percentage*100
...@@ -276,15 +218,9 @@ donut_d4_final <- (donut_d4_AAV9 | donut_d4_UT) ...@@ -276,15 +218,9 @@ donut_d4_final <- (donut_d4_AAV9 | donut_d4_UT)
## markers ## markers
for (t in tp) { for (t in tp) {
obj_t <- get(paste0(t,"_final")) obj_t <- get(paste0(t,"_final"))
neurons <- subset( neurons <- subset(obj_t, Cell_Type=="Neurons" | Cell_Type=="Immature.Neurons")
obj_t, astro <- subset(obj_t, Cell_Type=="Astrocytes" | Cell_Type=="Immature.Astrocytes")
RNA_snn_h.orig.ident_res.1.2=="Neurons" | RNA_snn_h.orig.ident_res.1.2=="Immature.Neurons") oligo <- subset(obj_t, Cell_Type=="Oligodendrocytes" | Cell_Type=="OPC")
astro <- subset(
obj_t,
RNA_snn_h.orig.ident_res.1.2=="Astrocytes" | RNA_snn_h.orig.ident_res.1.2=="Immature.Astrocytes")
oligo <- subset(
obj_t,
RNA_snn_h.orig.ident_res.1.2=="Oligodendrocytes" | RNA_snn_h.orig.ident_res.1.2=="OPC")
assign(paste0("neurons_",t,"_obj_x_markers"), neurons) assign(paste0("neurons_",t,"_obj_x_markers"), neurons)
assign(paste0("astro_",t,"_obj_x_markers"), astro) assign(paste0("astro_",t,"_obj_x_markers"), astro)
assign(paste0("oligo_",t,"_obj_x_markers"), oligo) assign(paste0("oligo_",t,"_obj_x_markers"), oligo)
......
scRNAseq_organoids/seurat_analysis.R
View file @ 31eabd95
...@@ -11,70 +11,21 @@ suppressPackageStartupMessages(library(ggrepel)) ...@@ -11,70 +11,21 @@ suppressPackageStartupMessages(library(ggrepel))
suppressPackageStartupMessages(library(clusterProfiler)) suppressPackageStartupMessages(library(clusterProfiler))
suppressPackageStartupMessages(library(RColorBrewer)) suppressPackageStartupMessages(library(RColorBrewer))
## load seurat data ## load seurat data
h48_final <- readRDS("data/organoids_h48_final.rds") tp <- c("h48","d5")
d5_final <- readRDS("data/organoids_D5_final.rds") for (t in tp) {
Idents(h48_final) <- "RNA_snn_h.orig.ident_res.0.6" sc_obj <- readRDS(paste0("data/organoids_",t,"_final.rds"))
Idents(d5_final) <- "RNA_snn_h.orig.ident_res.0.6" Idents(sc_obj) <- "Cell_Type"
h48_final@meta.data[['orig.ident']] <- recode_factor( sc_obj$Cell_type <- factor(
h48_final@meta.data[["orig.ident"]], sc_obj$Cell_Type,
"Sample1" = "AAV2", levels = c("Astrocytes",
"Sample3" = "AAV9", "Neurons",
"Sample5" = "Spk", "Oligodendrocytes",
"Sample7" = "UT") "Immature.Astrocytes",
h48_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]] <- recode_factor( "Immature.Neurons",
h48_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]], "OPC",
"0" = "Astrocytes", "Undifferentiated"))
"1" = "Immature.Astrocytes", assign(paste0(t,"_final"), sc_obj)
"2" = "Undifferentiated", }
"3" = "Neurons",
"4" = "Immature.Neurons",
"5" = "OPC",
"6" = "Astrocytes",
"7" = "Oligodendrocytes",
"8" = "Undifferentiated",
"9" = "Astrocytes",
"10" = "Undifferentiated",
"11" = "Undifferentiated")
h48_final$RNA_snn_h.orig.ident_res.0.6 <- factor(
h48_final$RNA_snn_h.orig.ident_res.0.6,
levels = c("Astrocytes",
"Neurons",
"Oligodendrocytes",
"Immature.Astrocytes",
"Immature.Neurons",
"OPC",
"Undifferentiated"))
d5_final@meta.data[['orig.ident']] <- recode_factor(
d5_final@meta.data[["orig.ident"]],
"Sample2" = "AAV2",
"Sample4" = "AAV9",
"Sample6" = "Spk",
"Sample8" = "UT")
d5_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]] <- recode_factor(
d5_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]],
"0" = "Astrocytes",
"1" = "Immature.Neurons",
"2" = "Undifferentiated",
"3" = "Neurons",
"4" = "Oligodendrocytes",
"5" = "Undifferentiated",
"6" = "Undifferentiated",
"7" = "Immature.Astrocytes",
"8" = "OPC",
"9" = "Immature.Neurons",
"10" = "Astrocytes",
"11" = "Undifferentiated",
"12" = "Astrocytes",
"13" = "Neurons")
d5_final$RNA_snn_h.orig.ident_res.0.6 <- factor(
d5_final$RNA_snn_h.orig.ident_res.0.6,
levels = c("Astrocytes",
"Neurons",
"Oligodendrocytes",
"Immature.Astrocytes",
"Immature.Neurons",
"OPC",
"Undifferentiated"))
ggplotColours <- function(n = 6, h = c(0, 360) + 15){ ggplotColours <- function(n = 6, h = c(0, 360) + 15){
if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65) hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
...@@ -82,7 +33,7 @@ ggplotColours <- function(n = 6, h = c(0, 360) + 15){ ...@@ -82,7 +33,7 @@ ggplotColours <- function(n = 6, h = c(0, 360) + 15){
## umap ## umap
DimPlot(object = h48_final, DimPlot(object = h48_final,
reduction = "umap.harmony.orig.ident", reduction = "umap.harmony.orig.ident",
group.by = "RNA_snn_h.orig.ident_res.0.6", group.by = "Cell_Type",
label = T, label = T,
repel = T, repel = T,
pt.size = 0.2, pt.size = 0.2,
...@@ -97,7 +48,7 @@ DimPlot(object = h48_final, ...@@ -97,7 +48,7 @@ DimPlot(object = h48_final,
values = ggplotColours(n = 7)) values = ggplotColours(n = 7))
## donut plot ## donut plot
ggplot_object_clusters_labels <- melt(table(h48_final$RNA_snn_h.orig.ident_res.0.6)) ggplot_object_clusters_labels <- melt(table(h48_final$Cell_Type))
colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count") colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count)) totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
totalB$percentage <- totalB$percentage*100 totalB$percentage <- totalB$percentage*100
...@@ -149,10 +100,10 @@ calc_helper <- function(object,genes){ ...@@ -149,10 +100,10 @@ calc_helper <- function(object,genes){
}else{return(NA)} }else{return(NA)}
} }
aav9 <- subset(h48_final, orig.ident=="AAV9") aav9 <- subset(d5_final, orig.ident=="AAV9")
gfp_percentage_cell_type <- PrctCellExpringGene(aav9, gfp_percentage_cell_type <- PrctCellExpringGene(aav9,
genes ="GFP", genes ="GFP",
group.by = "RNA_snn_h.orig.ident_res.0.6") group.by = "Cell_Type")
gfp_percentage_cell_type$Percentage <- gfp_percentage_cell_type$Cell_proportion*100 gfp_percentage_cell_type$Percentage <- gfp_percentage_cell_type$Cell_proportion*100
gfp_percentage_cell_type <- gfp_percentage_cell_type[,c(4,3)] gfp_percentage_cell_type <- gfp_percentage_cell_type[,c(4,3)]
colnames(gfp_percentage_cell_type) <- c("Percentage", "Cell_Types") colnames(gfp_percentage_cell_type) <- c("Percentage", "Cell_Types")
...@@ -190,49 +141,14 @@ ggplot(data=gfp_percentage_cell_type, ...@@ -190,49 +141,14 @@ ggplot(data=gfp_percentage_cell_type,
axis.text.y = element_text(size=18*96/72), axis.text.y = element_text(size=18*96/72),
legend.position = "none", legend.position = "none",
aspect.ratio = 1.1) aspect.ratio = 1.1)
h48_aav9_ut <- subset(h48_final, orig.ident == "AAV9" | orig.ident == "UT")
d4_aav9_ut <- subset(d4_final, orig.ident == "AAV9" | orig.ident == "UT")
DimPlot(object = h48_aav9_ut,
reduction = "umap.harmony.orig.ident",
group.by = "RNA_snn_h.orig.ident_res.1.2",
split.by = "orig.ident",
label = F,
repel = T,
pt.size = 0.2,
label.size = 6,
ncol = 2) +
scale_color_manual(labels=c("Astrocytes",
"Neurons",
"Oligodendrocytes",
"Immature Astrocytes",
"Immature Neurons",
"OPC",
"Undifferentiated"),
values = ggplotColours(n = 7))
DimPlot(object = d4_aav9_ut,
reduction = "umap.harmony.orig.ident",
group.by = "RNA_snn_h.orig.ident_res.1.2",
split.by = "orig.ident",
label = F,
repel = T,
pt.size = 0.2,
label.size = 6,
ncol = 2) +
scale_color_manual(labels=c("Astrocytes",
"Neurons",
"Oligodendrocytes",
"Immature Astrocytes",
"Immature Neurons",
"OPC",
"Undifferentiated"),
values = ggplotColours(n = 7))
samples <- c("AAV9", "UT") samples <- c("AAV9", "UT")
tp <- c("h48","d4") tp <- c("h48","d5")
for (t in tp) { for (t in tp) {
obj_t <- get(paste0(t,"_final")) obj_t <- get(paste0(t,"_final"))
for (i in samples) { for (i in samples) {
obj_i <- subset(obj_t, orig.ident == i) obj_i <- subset(obj_t, orig.ident == i)
ggplot_object_clusters_labels <- melt(table(obj_i$RNA_snn_h.orig.ident_res.1.2)) ggplot_object_clusters_labels <- melt(table(obj_i$Cell_Type))
colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count") colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count)) totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
totalB$percentage <- totalB$percentage*100 totalB$percentage <- totalB$percentage*100
...@@ -264,41 +180,88 @@ for (t in tp) { ...@@ -264,41 +180,88 @@ for (t in tp) {
assign(paste0("donut_",t,"_", i), donut_i) assign(paste0("donut_",t,"_", i), donut_i)
} }
} }
donut_h48_final <- (donut_h48_AAV9 | donut_h48_UT) donut_final <- (donut_h48_UT | donut_h48_AAV9 | donut_d5_UT | donut_d5_AAV9)
donut_d4_final <- (donut_d4_AAV9 | donut_d4_UT)
## markers ## markers
for (t in tp) { for (t in tp) {
obj_t <- get(paste0(t,"_final")) obj_t <- get(paste0(t,"_final"))
neurons <- subset( neurons <- subset(obj_t, Cell_Type=="Neurons" | Cell_Type=="Immature.Neurons")
obj_t, astro <- subset(obj_t, Cell_Type=="Astrocytes" | Cell_Type=="Immature.Astrocytes")
RNA_snn_h.orig.ident_res.1.2=="Neurons" | RNA_snn_h.orig.ident_res.1.2=="Immature.Neurons") oligo <- subset(obj_t, Cell_Type=="Oligodendrocytes" | Cell_Type=="OPC")
astro <- subset(
obj_t,
RNA_snn_h.orig.ident_res.1.2=="Astrocytes" | RNA_snn_h.orig.ident_res.1.2=="Immature.Astrocytes")
oligo <- subset(
obj_t,
RNA_snn_h.orig.ident_res.1.2=="Oligodendrocytes" | RNA_snn_h.orig.ident_res.1.2=="OPC")
assign(paste0("neurons_",t,"_obj_x_markers"), neurons) assign(paste0("neurons_",t,"_obj_x_markers"), neurons)
assign(paste0("astro_",t,"_obj_x_markers"), astro) assign(paste0("astro_",t,"_obj_x_markers"), astro)
assign(paste0("oligo_",t,"_obj_x_markers"), oligo) assign(paste0("oligo_",t,"_obj_x_markers"), oligo)
} }
for (k in ls()[grep("_obj_x_markers", ls())]) {
cells <- get(k) cell_types <- c("neurons_h48", "astro_h48", "oligo_h48",
degs <- FindMarkers(object = cells, "neurons_d5", "astro_d5", "oligo_d5")
ident.1 = "AAV9", samples_degs <- c("AAV9", "AAV2", "Spk")
ident.2 = "UT",
group.by = "orig.ident", for (k in cell_types) {
only.pos = FALSE, cells <- get(paste0(k,"_obj_x_markers"))
min.pct = 0, for (j in samples_degs) {
test.use = "wilcox", print(paste("Extracting differentially expressed genes of",j,"vs UT"))
logfc.threshold = 0, degs <- FindMarkers(cells,
min.cells.group = 0) ident.1 = j,
assign(paste0(k,"_aav9_vs_ut_degs"), degs) ident.2 = "UT",
group.by = "orig.ident",
only.pos = FALSE,
min.pct = 0.1,
test.use = "wilcox",
logfc.threshold = 0.1)
assign(paste(k,j,"vs_UT_degs",sep = "_"), degs, envir = .GlobalEnv)
degs_all <- FindMarkers(object = cells,
ident.1 = j,
ident.2 = "UT",
group.by = "orig.ident",
only.pos = FALSE,
min.pct = 0,
test.use = "wilcox",
logfc.threshold = 0,
min.cells.group = 0)
assign(paste(k,j,"vs_UT_degs_all",sep = "_"), degs_all, envir = .GlobalEnv)
## volcano plot
print(paste("Making volcano plot for DEGs of...", j, " vs UT" ))
logfc.thr <- 0
adj.pval.thr <- 0.05
data <- data.frame(gene = row.names(degs),
pval = -log10(degs$p_val_adj),
lfc = degs$avg_log2FC)
data <- na.omit(data)
data <- mutate(data,color=case_when(data$lfc>logfc.thr & data$pval> -log10(adj.pval.thr)~"Overexpressed",
data$lfc< -logfc.thr & data$pval> -log10(adj.pval.thr)~"Underexpressed",
abs(data$lfc)<logfc.thr|data$pval< -log10(adj.pval.thr)~"NonSignificant"))
data <- data[order(data$pval, decreasing = TRUE),]
highl_up <- head(subset(data, color != "NonSignificant" & lfc > 0), 10)
highl_down <- head(subset(data, color != "NonSignificant" & lfc < 0), 10)
highl <- rbind(highl_up, highl_down)
sig <- subset(degs, degs$p_val_adj<0.05)
n_sig <- nrow(sig)
vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
geom_point(size = 3, alpha = 0.8, na.rm = T) +
geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") +
geom_label_repel(data = highl,
aes(label = gene),
show.legend = FALSE,
size = 6,
family = "Arial") +
xlab(expression(log[2]("Fold Change"))) +
ylab(expression(-log[10]("adjusted p-value"))) +
scale_color_manual(name = "Expression",
values = c(Overexpressed = "red",
Underexpressed = "royalblue3",
NonSignificant = "black")) +
ggtitle(paste("Significant genes = ", n_sig)) +
theme_bw()
assign(paste(k,j,"vs_UT_volcano",sep = "_"), vol, envir = .GlobalEnv)
}
} }
rm(degs_all)
## gsea ## gsea
h_gmt <- read.gmt("data/h.all.v7.2.symbols.gmt") h_gmt <- read.gmt("data/h.all.v7.2.symbols.gmt")
for (i in ls()[grep("aav9_vs_ut_degs", ls())]) { for (g in ls()[grep("degs_all", ls())]) {
gene_res <- get(i) gene_res <- get(g)
logfc_symbol <- as.vector(gene_res$avg_log2FC) logfc_symbol <- as.vector(gene_res$avg_log2FC)
names(logfc_symbol) <- row.names(gene_res) names(logfc_symbol) <- row.names(gene_res)
logfc_symbol <- sort(logfc_symbol, decreasing = TRUE) logfc_symbol <- sort(logfc_symbol, decreasing = TRUE)
...@@ -307,60 +270,63 @@ for (i in ls()[grep("aav9_vs_ut_degs", ls())]) { ...@@ -307,60 +270,63 @@ for (i in ls()[grep("aav9_vs_ut_degs", ls())]) {
nPerm = 10000, nPerm = 10000,
pvalueCutoff = 1) pvalueCutoff = 1)
contr.gsea.symb <- as.data.frame(contr.gsea) contr.gsea.symb <- as.data.frame(contr.gsea)
assign(paste0("GSEA_",i), contr.gsea.symb) assign(paste0("GSEA_",g), contr.gsea.symb)
} }
## GSEA heatmaps ## GSEA heatmaps
type <- c("astro", "neurons", "oligo") type <- c("astro", "neurons", "oligo")
samples_degs <- c("AAV9", "AAV2", "Spk")
for (p in type) { for (p in type) {
hallmark.full <- data.frame() for (j in samples_degs) {
for (ds in ls()[grep(paste0("GSEA_", p), ls())]) { hallmark.full <- data.frame()
ds.t <- get(ds) for (ds in ls()[grep(paste0("GSEA_.*",p,".*",j), ls())]) {
ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalues"), drop = FALSE] ds.t <- get(ds)
ds.t.filt$Dataset <- ds ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalue"), drop = FALSE]
if (nrow(hallmark.full) == 0) { ds.t.filt$Dataset <- ds
hallmark.full <- ds.t.filt if (nrow(hallmark.full) == 0) {
} else { hallmark.full <- ds.t.filt
hallmark.full <- rbind(hallmark.full, ds.t.filt) } else {
hallmark.full <- rbind(hallmark.full, ds.t.filt)
}
} }
hallmark.full.filt <- hallmark.full[, c("ID", "NES", "Dataset", "qvalue")]
hallmark.full.filt$ID <- gsub(x = hallmark.full.filt$ID, "HALLMARK_", "")
hallmark.full.filt$sig <- ifelse(
hallmark.full.filt$qvalue <= 0.05 & hallmark.full.filt$qvalue > 0.01, '*',
ifelse(hallmark.full.filt$qvalue <= 0.01 & hallmark.full.filt$qvalue > 0.001, '**',
ifelse(hallmark.full.filt$qvalue <= 0.001 & hallmark.full.filt$qvalue > 0.0001, '***',
ifelse(hallmark.full.filt$qvalue <= 0.0001, '****', ''))))
hallmark.order <- hallmark.full.filt %>%
group_by(ID) %>%
summarise(Pos = sum(NES))
hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
hallmark.full.filt.tt <- reshape2::melt(reshape2::dcast(hallmark.full.filt,
ID ~ Dataset,
value.var="NES"), id.vars = c("ID"))
colnames(hallmark.full.filt.tt) <- c("ID", "Dataset", "NES")
hallmark.full.filt.tt$ID <- factor(hallmark.full.filt.tt$ID, levels = hallmark.order.terms$ID)
hallmark.full.filt.tt <- merge(hallmark.full.filt.tt,
hallmark.full.filt[, c("ID", "Dataset", "sig")],
by = c("ID", "Dataset"),
all.x = TRUE)
labels <- c("Day4", "Day2")
p.hallmark <- ggplot(hallmark.full.filt.tt,
aes(x = Dataset, y = ID)) +
geom_tile(aes(fill = NES), colour = "white") +
geom_text(aes(label = paste(sig)), size=4*96/72, fontface = "bold") +
scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
limits = c(-3.5, 3.5),
na.value = "white") +
ylab("") +
xlab("") +
coord_fixed(ratio = 0.6) +
scale_x_discrete(labels = labels) +
theme_bw(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size = 18*96/72),
axis.text.y = element_text(face = "bold", size = 12*96/72),
legend.text = element_text(face = "bold", size = 14*96/72),
legend.title = element_text(face = "bold", size = 12*96/72))
assign(paste("p.hallmark", p,j,sep = "_"), p.hallmark)
} }
hallmark.full.filt <- hallmark.full[, c("ID", "NES", "Dataset", "qvalues")]
hallmark.full.filt$ID <- gsub(x = hallmark.full.filt$ID, "HALLMARK_", "")
hallmark.full.filt$sig <- ifelse(
hallmark.full.filt$qvalues <= 0.05 & hallmark.full.filt$qvalues > 0.01, '*',
ifelse(hallmark.full.filt$qvalues <= 0.01 & hallmark.full.filt$qvalues > 0.001, '**',
ifelse(hallmark.full.filt$qvalues <= 0.001 & hallmark.full.filt$qvalues > 0.0001, '***',
ifelse(hallmark.full.filt$qvalues <= 0.0001, '****', ''))))
hallmark.order <- hallmark.full.filt %>%
group_by(ID) %>%
summarise(Pos = sum(NES))
hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
hallmark.full.filt.tt <- reshape2::melt(reshape2::dcast(hallmark.full.filt,
ID ~ Dataset,
value.var="NES"), id.vars = c("ID"))
colnames(hallmark.full.filt.tt) <- c("ID", "Dataset", "NES")
hallmark.full.filt.tt$ID <- factor(hallmark.full.filt.tt$ID, levels = hallmark.order.terms$ID)
hallmark.full.filt.tt <- merge(hallmark.full.filt.tt,
hallmark.full.filt[, c("ID", "Dataset", "sig")],
by = c("ID", "Dataset"),
all.x = TRUE)
labels <- c("Day4", "Day2")
p.hallmark <- ggplot(hallmark.full.filt.tt,
aes(x = Dataset, y = ID)) +
geom_tile(aes(fill = NES), colour = "white") +
geom_text(aes(label = paste(sig)), size=4*96/72, fontface = "bold") +
scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
limits = c(-3.5, 3.5),
na.value = "white") +
ylab("") +
xlab("") +
coord_fixed(ratio = 0.6) +
scale_x_discrete(labels = labels) +
theme_bw(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size = 18*96/72),
axis.text.y = element_text(face = "bold", size = 12*96/72),
legend.text = element_text(face = "bold", size = 14*96/72),
legend.title = element_text(face = "bold", size = 12*96/72))
assign(paste0("p.hallmark_", p), p.hallmark)
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment