基因富集分析

immu_data<-load('E:\\dt.Tcell.Rdata')

library(Seurat)
immu_data<-dt.Tcell

table(Idents(immu_data))
q

expr_all_ <- immu_data
Idents(expr_all_) <- expr_all$sample.ident
expr_all <- expr_all_

library(clusterProfiler)
library(org.Mm.eg.db)
#gl <- hallmarks@.Data[[26]]@geneIds
#rowgene_h <- readRDS("D:/Gu_lab/Liver/Data/mm2hm.RDS")
#gl <- rowgene_h$MGI.symbol[which(rowgene_h$HGNC.symbol%in%gl)]

gl <- intersect(ind,rownames(expr_all))

p1 <- VlnPlot(expr_all_,gl,pt.size = 0, cols = cols2)
ggsave(filename = paste0(merge_savePath,"Kegg_mTor_gene_source.png"), p1,
width = 10, height = 28, dpi = 900,limitsize = F)

p1 <- VlnPlot(expr_all,gl,pt.size = 0,cols = cols1)

ggsave(filename = paste0(merge_savePath,"Kegg_mTor_gene.png"), p1,
width = 15, height = 28, dpi = 500,limitsize = F)

p1 <- FeaturePlot(expr_all,gl)

ggsave(filename = paste0(merge_savePath,"Other_gene.png"), p1,
width = 19, height = 200/13, dpi = 500,limitsize = F)

expr_all_[['lotus']] <- colMeans(expr_all_@assays$RNA@data[gl,])
cols = as.array(getDefaultColors(length(table(Idents(expr_all_)))))
rownames(cols) = names(table(Idents(expr_all_)))
p2 <- VlnPlot(expr_all_,'lotus',pt.size = 0) +
scale_fill_manual(values = cols,
guide = guide_legend(override.aes = list(size = 3),
keywidth = 0.1,
keyheight = 0.15,
default.unit = "inch"))
ggsave(filename = paste0(save_Path,"Proliferating_source.png"), p2,
width = 5, height = 4, dpi = 300)

##GSVA

library(GSVA)
library(GSEABase)
library(msigdbr)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(limma)

pw <- unique(na.omit(gl))

pathway_list <- list(pw)

names(pathway_list) <- c("lotus")
rm(expr_all_)
gsva_matrix_BD <- gsva(as.matrix(expr_all@assays$RNA@scale.data), pathway_list,method='gsva')
write.csv(gsva_matrix_BD,file = "gsva_matrix_BD.csv")

expr_all[['lotus_gsva']] <- gsva_matrix_BD[1,]
cols = as.array(getDefaultColors(length(table(Idents(expr_all)))))
rownames(cols) = names(table(Idents(expr_all)))
p2 <- VlnPlot(expr_all,'lotus_gsva',pt.size = 0) +
scale_fill_manual(values = cols,
guide = guide_legend(override.aes = list(size = 3),
keywidth = 0.1,
keyheight = 0.15,
default.unit = "inch"))
ggsave(filename = paste0(save_Path,"Proliferating_source.png"), p2,
width = 5, height = 4, dpi = 300)

#GSEA
library(clusterProfiler)
markers <- FindAllMarkers(expr_all)
markers <- markers[which(markers$avg_log2FC>0.05),]
id <- names(table(markers$cluster))
score <- c()
count <- c()
pw <- data.frame(cellName = "lotus", geneID = sample(gl,500))
pw <- as_tibble(pw)
for(i in id){
gl_ <- markers[which(markers$cluster==i),]
gl_ <- gl_[sort(gl_$avg_log2FC,index.return=TRUE)$ix,]
gsea_ <- enricher(gl_$gene,TERM2GENE = pw,pvalueCutoff = 1)
score <- c(score,gsea_@result$pvalue)
count <- c(count,gsea_@result$Count)
}
names(count) <- id
count

###################function##################
#https://yulab-smu.top/biomedical-knowledge-mining-book/universal-api.html