Commit 5ffb27de authored by Tesset1's avatar Tesset1
Browse files

add script

parent 39bc00b5
########################
## Single cell analysis
########################
library(Seurat)
library(clusterProfiler)
library(patchwork)
library(ggplot2)
library(RColorBrewer)
library(scales)
library(grid)
library(gridExtra)
## full data all the condition
readRDS('full_csc_2023.rds')
## senescence lists
gmt.obj <- read.gmt('senescence_list_paper.gmt')
gmt.obj$term <- factor(gmt.obj$term)
table(gmt.obj$term)
ll <- unique(gmt.obj$term)
full_obj$CSC_clusters <- ifelse(full_obj$CSC_celltype == 'HSC_0', 'HSC-enriched_1',
ifelse(full_obj$CSC_celltype == 'HSC_1', 'HSC-enriched_1',
ifelse(full_obj$CSC_celltype == 'HSC_2', 'HSC-enriched_2',
ifelse(full_obj$CSC_celltype == 'HSC_3', 'HSC-enriched_1',
ifelse(full_obj$CSC_celltype == '1', 'Myeloid-biased_prog',
ifelse(full_obj$CSC_celltype == '3', 'Erythroid-biased_prog',
ifelse(full_obj$CSC_celltype == '4', 'Neutrophil-biased_prog',
ifelse(full_obj$CSC_celltype == '5', 'PreB/Mono/DC_prog','unknown'
))))))))
full_obj$CSC_clusters <- factor(full_obj$CSC_clusters, levels=c('HSC-enriched_1', 'HSC-enriched_2',
'PreB/Mono/DC_prog', 'Myeloid-biased_prog',
'Neutrophil-biased_prog','Erythroid-biased_prog'
))
table( colnames(mix3_obj_new) == colnames(full_obj))
full_obj$neworignames <- ifelse(full_obj$orig.ident == 'RNPneg', 'RNP NEG',
ifelse(full_obj$orig.ident == 'GTGFPneg', 'GFP-',
ifelse(full_obj$orig.ident == 'GTGFPpos', 'GFP+',
ifelse(full_obj$orig.ident == 'RNPhs', 'RNP HS',full_obj$orig.ident)))))
Idents(full_obj) <- 'neworignames'
full_obj <- subset(full_obj, idents = c("GFP+", "GFP-","RNP NEG", "RNP HS"))
full_obj$neworignames <- factor(full_obj$neworignames, levels =c("RNP NEG","RNP HS",
"GFP-","GFP+"))
Idents(full_obj) <- 'CSC_clusters'
for (l in 1:length(ll)) {
group <-ll[l]
nn <- gsub('\\)','',gsub('\\(','',gsub(' ','_',group)))
nn <- gsub('-','_',nn)
glist <- gmt.obj[gmt.obj$term == group,'gene']
gmt_features <- list(glist)
full_obj <- AddModuleScore(
object = full_obj,
features = gmt_features, ##
name = nn,
)}
clcol <- c("RNP NEG"="black", "RNP HS"="gray36", "GFP-"="grey", "GFP+"="green3")
Idents(full_obj) <- 'CSC_clusters'
p <- list()
sasp <- list()
for (n in levels(full_obj$CSC_clusters)){
sub <- subset(full_obj, idents = n)
p[[n]] <- VlnPlot(object = sub, features = c("FRIDMAN_SENESCENCE_UP1"),
group.by ="neworignames", col = clcol, pt.size = 0,sort = F)+
geom_boxplot(width=0.2,outlier.shape = NA)+ theme(legend.position="none",
axis.title.x = element_blank())+
ggtitle(paste0(gsub('_', ' ',n)))
sasp[[n]] <- VlnPlot(object = sub, features = c("SASP_SCHLEICH1"),
group.by ="neworignames", col = clcol, pt.size = 0,sort = F)+
geom_boxplot(width=0.2,outlier.shape = NA)+ theme(legend.position="none",
axis.title.x = element_blank())+
ggtitle(paste0(gsub('_', ' ',n)))
}
## get legend just from one
lg_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
mylegend<- lg_legend(VlnPlot(object = sub, features = c("FRIDMAN_SENESCENCE_UP1"),
group.by ="neworignames", col = clcol, pt.size = 0,sort = F)+
geom_boxplot(width=0.2,outlier.shape = NA))
svg("senescence_violin.svg", width =20, height = 4, pointsize=12)
grid.arrange(p$`HSC-enriched_1`, p$`HSC-enriched_2`, p$`PreB/Mono/DC_prog`,
p$`Myeloid-biased_prog`, p$`Neutrophil-biased_prog`, p$`Erythroid-biased_prog`, mylegend,
nrow = 1,top = textGrob("FRIDMAN SENESCENCE UP",gp=gpar(fontsize=20,font=3)))
dev.off()
svg("sasp_violin.svg", width =20, height = 4, pointsize=12)
grid.arrange(sasp$`HSC-enriched_1`, sasp$`HSC-enriched_2`, sasp$`PreB/Mono/DC_prog`,
sasp$`Myeloid-biased_prog`, sasp$`Neutrophil-biased_prog`, sasp$`Erythroid-biased_prog`, mylegend,
nrow = 1,top = textGrob("SASP SCHLEICH",gp=gpar(fontsize=20,font=3)))
dev.off()
## RIDGE plot
p <- RidgePlot(object = full_obj, features = c("FRIDMAN_SENESCENCE_UP1"),group.by = "neworignames", cols = clcol)+ theme(legend.position="none",
axis.title.x = element_blank(),
axis.title.y = element_blank())
svg("senescence_ridge.svg", width =11, height = 6, pointsize=12)
p + ggtitle("FRIDMAN SENESCENCE UP")
dev.off()
## final sasp
p <- RidgePlot(object = full_obj, features = c("SASP_SCHLEICH1"),group.by = "neworignames", cols = clcol)+ theme(legend.position="none",
axis.title.x = element_blank(),
axis.title.y = element_blank())
svg("sasp_SCHLEICH_ridge.svg", width =11, height = 6, pointsize=12)
p + ggtitle("SASP SCHLEICH")
dev.off()
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