GSEABase | clusterProfiler | GSEA富集分析 | MSigDB
2023年09月01日
log2FC <- log2(rowMeans(TPM[,high.samples])+1) - log2(rowMeans(TPM[,low.samples])+1) geneList <- log2FC geneList <- sort(geneList, decreasing = T) egmt <- clusterProfiler::GSEA(geneList[geneList!=0], TERM2GENE = geneset, minGSSize = 1, pvalueCutoff = 1, nPermSimple = 10000, verbose = F) options(repr.plot.width=8, repr.plot.height=4) enrichplot::gseaplot2(egmt, geneSetID = "Enterocyte-marker genes (Haber et al., 2017)", pvalue_table = T, subplots = 1:2, ES_geom="line", color = "red")
term | gene |
gene set name | gene name |
2023年08月22日
美图参考:http://localhost:17435/notebooks/data_center/DB/DB-TCGA-CCLE-GTEx.ipynb
2022年10月07日
用MSigDB的数据集来做GSEA分析,这是常规可靠的分析,能上CNS。
最常见的第一种hallmark数据集
第二种就是GOBP
参考:Data_center/analysis/ApcKO_multiomics/3.DEG_pathways.ipynb
注意数据集需要clean,基因个数不能少于5,否则没有富集得分,程序会报错。
核心的函数是
tmp.NES <- clusterProfiler::GSEA(geneList2, TERM2GENE = gmt.list[[j]], minGSSize = 1, pvalueCutoff = pvalueCutoff, verbose = F)
MSigDB可以用msigdbr来读取
h.gmt <- msigdbr(species = species, category = "H") %>% dplyr::select(gs_name, gene_symbol) # entrez_gene
msigdbr_collections()
参考:
2022年09月07日
太厉害了,可以做自定义的gene set的GSEA分析。
参考:
- 处理gmt文件的一些技巧
- /project/Data_center/analysis/Perturb_seq_Sethi/cellranger/DEG_pathway.ipynb
结果非常好,是可以发CNS的水平。
2022年08月31日
参考:Data_center/analysis/Perturb_seq_Sethi/cellranger/DEG_pathway.ipynb
其实我很早之前就写了个gsea的函数,也是用clusterProfiler来做,但是貌似结果都不是很理想,老板也不爱用这个分析,所以我也就一直没用。
今天看着别人的代码,把自己的函数又重新调整了一下。
参考:生信技能树 - 代码有所更新
获取单细胞测试数据
# devtools::install_github("satijalab/seurat-data") library(SeuratData) # AvailableData() # InstallData("pbmc3k.SeuratData") data(pbmc3k) exp <- pbmc3k@assays$RNA@data dim(exp) # exp[1:5,1:5] table(is.na(pbmc3k$seurat_annotations)) table(pbmc3k$seurat_annotations) library(Seurat) pbmc3k@active.ident <- pbmc3k$seurat_annotations table(pbmc3k@active.ident) deg <- FindMarkers(pbmc3k, ident.1 = "Naive CD4 T", ident.2 = "B") # head(deg) dim(deg)
GSEA分析
# if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("GSEABase") library(GSEABase) library(ggplot2) library(clusterProfiler) library(org.Hs.eg.db) # API 1 geneList <- deg$avg_logFC names(geneList) <- toupper(rownames(deg)) geneList <- sort(geneList, decreasing = T) head(geneList) # API 2 # gmtfile <- "../EllyLab//human/singleCell/MsigDB/msigdb.v7.4.symbols.gmt" gmtfile <- "../EllyLab//human/singleCell/MsigDB/c5.go.bp.v7.4.symbols.gmt" geneset <- read.gmt(gmtfile) length(unique(geneset$ont)) egmt <- GSEA(geneList, TERM2GENE = geneset, minGSSize = 1, pvalueCutoff = 0.99, verbose = F) # head(egmt) gsea.out.df <- egmt@result # gsea.out.df # head(gsea.out.df$ID) library(enrichplot)
出图 - 基因数据足够才够漂亮
options(repr.plot.width=6, repr.plot.height=4) gseaplot2(egmt, geneSetID = "GOBP_ANTIGEN_RECEPTOR_MEDIATED_SIGNALING_PATHWAY", pvalue_table = T)
原始的图不够漂亮,优化可以参考阿汤哥的paper
这个图的代码不错,不知道他们paper里有没有分享。
2019 - Single-cell RNA-seq analysis reveals the progression of human osteoarthritis