富集分析-4分组代码

有分组，也就是两组对于某一个亚群的比较

参考视频：单细胞绘图之GSEA & GSVA_哔哩哔哩_bilibili

###11、GSVA分析  https://www.bilibili.com/video/BV1SD4y1q7zq/?spm_id_from=333.337.search-card.all.click&vd_source=7f9d77ce78399f44a43e8653e1c05c29
library(KEGGREST, quietly = TRUE)
library(tidyverse, quietly = TRUE)
BiocManager::install("GSVA")

library(clusterProfiler)
library(tidyverse)
library(ggplot2)
library(ggsci)
#BiocManager::install("org.Hs.eg.db") 
library(org.Hs.eg.db)



table(integrated@meta.data$seurat_clusters)

# 提取细胞的胞名
cells_NK2 <- subset(integrated@meta.data, seurat_clusters %in% c("4"))

# 把NK细胞从总数据集中提取出来
scRNA_sub = subset (integrated, cells = row.names(cells_NK2))  

view = scRNA_sub@meta.data


# 差异表达分析
sub.marker <- FindMarkers(integrated,group.by ='orig.ident', 
                            ident.1 = 'Use',ident.2 = 'Withdrawal',logfc.threshold = 0.01) 
#GSEA尽可能需要全部基因的logFC数据，所以logfc.threshold = 0.01，设置的小一些

sub.marker.sig = subset(sub.marker, p_val_adj<0.05 & avg_log2FC >0.5)


# GSEA

mydata = data.frame(Gene = row.names(sub.marker), logFC = sub.marker$avg_log2FC) %>%
  arrange(desc(logFC))  # 基因根据logFC排序

#把基因symbol转换成基因ID
genelist = bitr(mydata$Gene,fromType="SYMBOL",toType="ENTREZID", OrgDb="org.Hs.eg.db")%>%
                inner_join(mydata,by=c('SYMBOL'='Gene'))%>%
                arrange(desc(logFC))

               
#准备GSEA需要的输入文件
gsea_input = genelist$logFC
names(gsea_input) = genelist$ENTREZID
  
geneset = read.gmt( "D:\\atrial fibrillation 2022.11.03\\c7.all.v2023.1.Hs.entrez.gmt" )             

gg = GSEA(gsea_input, TERM2GENE = geneset, verbose = F,
          pvalueCutoff = 0.1,pAdjustMethod = "BH")


sortgg = gg[order(gg$NES,decreasing = T),]
sortgg = sortgg[sortgg$p.adjust < 0.05,]
save(gg,file = "NK_GSEA.rda")
load( "NK_GSEA.rda")


#根据实验设计、或者是P值的显著性、富集评分的显著性进行筛选
library(enrichplot)
geneset_plot = c("GSE22886_NAIVE_CD8_TCELL_VS_MONOCYTE_DN",
                 "GSE10325_LUPUS_CD4_TCELL_VS_LUPUS_MYELOID_DN",
                 "GSE22886_NAIVE_BCELL_VS_MONOCYTE_DN")
mycol = pal_nejm()(8)  # 图调色的参数
gseaplot2(gg,
          geneSetID = geneset_plot,
          color = mycol[c(1:3)],
          title = 'immunologic signature gene sets',
          rel_heights = c(1.3,0.3,0.6),
          pvalue_table = F)