大家好，欢迎来到IT知识分享网。
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()
参考：
- 如何使用clusterProfiler对MSigDB数据库进行富集分析
2022年09月07日
太厉害了，可以做自定义的gene set的GSEA分析。
参考：
- 处理gmt文件的一些技巧
- /project/Data_center/analysis/Perturb_seq_Sethi/cellranger/DEG_pathway.ipynb
结果非常好，是可以发CNS的水平。
2022年08月31日
参考：手把手教你用R做GSEA分析
参考：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
免责声明：本站所有文章内容,图片，视频等均是来源于用户投稿和互联网及文摘转载整编而成，不代表本站观点，不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益，请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权，请联系删除，如若转载，请注明出处：https://haidsoft.com/29103.html
