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