Step 7：Single-cell sequencing analysis and GSEA analysis

library(Seurat)

####################
logFCfilter=1               
adjPvalFilter=0.05   
ct=data.table::fread("GSE137665_exp_matrix.all.tsv.gz",data.table = F)
rownames(ct) <- ct$V1
ct <- ct[,-1]
data <- data.frame(colnames(ct))
data <- cbind(data, do.call(rbind, strsplit(as.character(data$colnames.ct.), split = "[-]")))

sce = CreateSeuratObject(counts =  ct ,
                         project =  "brain" ,
                         min.cells = 5,
                         min.features = 300,)

##############
sce.all <- sce
print(sce.all)
sce.all[["percent.mt"]] <- PercentageFeatureSet(object = sce.all, pattern = "^mt-")
sce.all=subset(x = sce.all, subset =nFeature_RNA > 300 &nFeature_RNA < 2000& percent.mt < 10)

library(ggplot2)
library(patchwork) 

# 设定颜色 (可选)
my_cols <- c("#4DBBD5", "#00A087") 

# 绘制
plot_list <- VlnPlot(object = sce.all,
                     features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
                     ncol = 3,
                     pt.size = 0,
                     cols = my_cols, 
                     combine = FALSE) 


plot_list <- lapply(plot_list, function(x) {
  x +
    geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
    theme_bw() +
    theme(
      axis.title.x = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1, color="black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      legend.position = "none"
    )
})


wrap_plots(plot_list, ncol = 3)
dev.off()


sce.all <- NormalizeData(sce.all)
sce.all <- FindVariableFeatures(sce.all)
sce.all <- ScaleData(sce.all)
sce.all <- RunPCA(sce.all)
sce.all <- FindNeighbors(sce.all, dims = 1:20, reduction = "pca")
sce.all <- FindClusters(sce.all, resolution = 0.5, cluster.name = "unintegrated_clusters")
sce.all <- RunUMAP(sce.all, dims = 1:20, reduction = "pca")

pdf(file = "UMAP.pdf",width=6,height=4.5)
DimPlot(sce.all)
dev.off()


###################################04.SingleR##################################################
new.cluster.ids <- c(
  "0" = "Endothelial",
  "1" = "Microglia",
  "2" = "Neuron Glutamatergic",
  "3" = "Neuron GABAergic",
  "4" = "Neuron GABAergic",
  "5" = "EPC",
  "6" = "Astrocyte",
  "7" = "Neuron Glutamatergic",
  "8" = "Astrocyte",
  "9" = "VSMCA",
  "10" = "Neuron GABAergic",
  "11" = "Neuron GABAergic",
  "12" = "Macrophage",
  "13" = "Neuron Other",
  "14" = "Neuron Other",
  "15" = "Endothelial",
  "16" = "Oligo",
  "17" = "Pericyte",
  "18" = "Astrocyte",
  "19" = "Astrocyte",
  "20" = "Microglia",
  "21" = "HbVC"
)
sce.all <- RenameIdents(sce.all, new.cluster.ids)
sce.all$cell_type <- Idents(sce.all)


UMAPPlot(sce.all,label = T)
genes <- list("Neuron" = c("Snap25","Dlk1","Chga"),
              "Astrocyte"=c( "Aldoc"),
              "EPC"=c("Tmem212","Rarres2"),
              "VSMCA"=c("Acta2"),
              "Pericyte"=c("Vtn"),
              "Endothelial"=c("Flt1"),
              "HbVC" = c("Alas2"),
              "Oligo" = c("Cldn11"),
              "Microglia"=c( "Cx3cr1"),
              "Macrophage"=c("Lyz2")
)
flat_genes <- unique(unlist(genes))

DotPlot1 <- DotPlot(sce.all, features = flat_genes, cols = c("lightgrey", "red")) + # 经典的灰红配色
  RotatedAxis() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, size = 1), # 加个黑框
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10))
ggplot2::ggsave("marker.pdf", plot = DotPlot1, device = "pdf",width=9,height=4)

Idents(sce.all) <- sce.all$cell_type

plot <- DimPlot(sce.all)
plot
ggplot2::ggsave("umap_celltype.pdf", plot = plot, device = "pdf",width=6.5,height=4.5)

save(sce.all,file="sce.all_celltype.RData")

############
hubgene <- c("Cdkn1a","Hspa5","Nr4a1")

DotPlot(object = sce.all, features = hubgene)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10))
##############
library(dplyr)
library(tidyr)
library(ggplot2)
library(pROC)
library(rms)

data_to_plot <- FetchData(sce.all, vars = c(hubgene, "cell_type"))
data_long <- data_to_plot %>%
  pivot_longer(cols = hubgene, names_to = "Gene", values_to = "Expression")
plot_data <- data_long %>%
  group_by(cell_type, Gene) %>%
  summarise(
    AvgExp = mean(expm1(Expression)), 
    PctExp = sum(Expression > 0) / n() * 100 
  ungroup()


plot_data <- plot_data %>%
  group_by(Gene) %>%
  mutate(AvgExpScaled = scale(AvgExp))

P <- ggplot(plot_data, aes(x = Gene, y = cell_type)) +
  geom_point(aes(size = PctExp, color = AvgExpScaled)) +
  scale_color_gradientn(colours = c("blue","lightgrey",  "red")) +
  theme_bw() +
  labs(size = "% Expressed", color = "Avg Exp (Scaled)", y = "Cluster", x = "Gene") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot2::ggsave("group_exp.pdf",plot = P, device = "pdf",width=5.5,height=3.5)
####################
library(dplyr)
library(ClusterGVis)
library(org.Hs.eg.db)
cell1 <- JoinLayers(sce.all)
save(cell1,file="cell1.RData")
sce.markers.all <- FindAllMarkers(cell1, logfc.threshold = 0.25, min.pct = 0.1)
save(sce.markers.all,file="sce.markers.all.RData")
sce.markers <- sce.markers.all%>%
  dplyr::group_by(cluster)%>%
  dplyr::top_n(n=30,wt=avg_log2FC)

common_genes <- intersect(rownames(cell1), sce.markers$gene)

st.data <- prepareDataFromscRNA(object=cell1,
                                diffData=sce.markers,
                                showAverage=TRUE)

enrich <- enrichCluster(object=st.data,
                        OrgDb= org.Mm.eg.db,
                        type="BP",
                        organism="hsa",
                        pvalueCutoff=0.5,
                        topn=5,
                        seed=123456)###BP,CC,MF
markGenes = unique(sce.markers$gene)[sample(1:length(unique(sce.markers$gene)),30,
                                            replace=F)]

pdf('term2.pdf',height = 8,width = 13,onefile = F)
visCluster(object = st.data,
           plot.type = "both",
           column_names_rot = 45,
           markGenes = markGenes,
           markGenes.side = "left",
           genes.gp = c('italic',fontsize = 12,col = "orange"),
           annoTerm.data = enrich)
dev.off()
library(org.Hs.eg.db)
library(stringr)
library(BiocGenerics)

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("clusterProfiler")
library(clusterProfiler)
library(enrichplot)
library(future)
library(future.apply)


setwd("-")
load("SD_exp_combined-human.Rdata")

exprSet=exp2


batch_cor <- function(gene){
  y <- as.numeric(exprSet[gene,])
  rownames <- rownames(exprSet)
  do.call(rbind,future_lapply(rownames, function(x){
    dd  <- cor.test(as.numeric(exprSet[x,]),y,type="spearman")
    data.frame(gene=gene,mRNAs=x,cor=dd$estimate,p.value=dd$p.value )
  }))
}


dd <- batch_cor("NR4A1")
gene <- dd$mRNAs
）
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
gene_df <- data.frame(cor=dd$cor,SYMBOL = dd$mRNAs)

gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

geneList <- gene_df$cor
names(geneList)=gene_df$SYMBOL

geneList=sort(geneList,decreasing = T)

dotplot_internal <- function(object, x = "GeneRatio", color = "pvalue",
                             showCategory=10, size=NULL, split = NULL,
                             font.size=12, title = "", orderBy="x", decreasing=TRUE) {
  
  colorBy <- match.arg(color, c("pvalue", "p.adjust", "qvalue"))
  
  if (x == "geneRatio" || x == "GeneRatio") {
    x <- "GeneRatio"
    if (is.null(size))
      size <- "Count"
  } else if (x == "count" || x == "Count") {
    x <- "Count"
    if (is.null(size))
      size <- "GeneRatio"
  } else if (is(x, "formula")) {
    x <- as.character(x)[2]
    if (is.null(size))
      size <- "Count"
  } else {
    if (is.null(size))
      size  <- "Count"
  }
  df <- fortify(object, showCategory = showCategory, split=split)
  
  if (orderBy !=  'x' && !orderBy %in% colnames(df)) {
    message('wrong orderBy parameter; set to default `orderBy = "x"`')
    orderBy <- "x"
  }
  
  if (orderBy == "x") {
    df <- dplyr::mutate(df, x = eval(parse(text=x)))
  }
  
  idx <- order(df[[orderBy]], decreasing = decreasing)
  df$Description <- factor(df$Description, levels=rev(unique(df$Description[idx])))
  ggplot(df, aes_string(x=x, y="Description", size=size, color=colorBy)) +
    geom_point() +
    scale_color_continuous(low="darkorange3", high="darkseagreen4", name = color, guide=guide_colorbar(reverse=TRUE)) +
    ylab(NULL) + ggtitle(title) + DOSE::theme_dose(font.size) + scale_size(range=c(3, 8))
  
}





GO<-read.gmt("c5.go.v2024.1.Hs.symbols.gmt")

GO <-GSEA(geneList,TERM2GENE = GO)

library(ggplot2)
library(enrichplot)


write.csv(data.frame(GO), file = "PRKAB1.csv", row.names = FALSE)


p1=gseaplot2(GO,
             geneSetID = 5,
             title=GO$Description[5],
             color="red", 
             base_size = 20, 
             subplots = 1:3, 
             pvalue_table = T) 
p1
ggsave("PRKAB1.pdf",width = 12,height = 10)

p2=gseaplot2(GO,
             geneSetID = 1:3,
             subplots = 1:3,
             ES_geom='line'
)
p2
ggsave("PRKAB1/p2.pdf",width = 12,height = 10)
num=5
p3=gseaplot2(GO, geneSetID = rownames(GO@result)[head(order(GO@result$enrichmentScore),num)])
p3
ggsave("PRKAB1/p3.pdf",width = 12,height = 10)
p4=gseaplot2(GO, geneSetID = rownames(GO@result)[tail(order(GO@result$enrichmentScore),num)])
p4
ggsave("PRKAB1/p4.pdf",width = 12,height = 10)
num=5
p5=gseaplot2(GO, geneSetID = rownames(GO@result)[c(head(order(GO@result$enrichmentScore),num),tail(order(GO@result$enrichmentScore),num))])
p5
ggsave("PRKAB1结果/p5.pdf",width = 12,height = 12)



#气泡图
p6=dotplot(GO,color="pvalue")
p6
ggsave("PRKAB1/p6.pdf",width = 10,height = 10)

#分类气泡图

p7=dotplot(GO,split=".sign")+facet_grid(~.sign)
p7
ggsave("PRKAB1/p7.pdf",width = 10,height = 10)