ixxmu
commented
1 year ago 奔向自由的富集分析 by 生信星球
今天是生信星球陪你的第919天
0.单细胞的marker基因
啥基因都可以，富集分析只需要一组基因名字。
不想要线粒体核糖体红细胞的基因，直接去掉。
rm(list = ls())
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
load("makers.Rdata")
k = grepl("^RP[SL]|^MT-|^HB[^(P)]",rownames(markers));table(k)
## k
## FALSE TRUE
## 16827 577
markers = markers[!k,]
g <- markers[markers$cluster == 6 &
markers$avg_log2FC > 3, 'gene' ]
g
## [1] "CTSW" "XCL2" "GZMK" "CD247" "CD2" "LCK"
## [7] "CD3E" "CST7" "PTPRCAP" "PRF1" "GZMA" "CXCR3"
## [13] "GZMH" "XCL1" "NKG7" "KLRD1" "ITM2A" "GZMM"
## [19] "KLRB1" "LAT" "CD3D" "ZAP70" "MATK" "CD96"
## [25] "KLRC1" "TBC1D10C" "CD7" "TMIGD2" "IL2RG" "HOPX"
## [31] "EVL" "IL2RB" "GIMAP7" "TRAF3IP3" "SELL" "CCL5"
## [37] "SEPT1" "CCR7" "ACAP1" "CD69" "CORO1A" "LIMD2"
## [43] "RHOH" "GZMB" "GNLY" "IL32" "PLAC8" "STK17A"
## [49] "S1PR4" "RAC2" "TNFRSF18" "ICAM3" "SEPT6" "CXCR4"
## [55] "PTPN7" "APOBEC3G" "PTPRC" "CDKN2D" "HLA-F" "HMGB2"
## [61] "PCSK7" "FDFT1" "SYNCRIP"
1.常规的富集分析
library(msigdbr)
ego <- enrichGO(
gene = g,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = "ALL"
)
dotplot(ego)
ego@result$Description[1:10]
## [1] "T cell activation"
## [2] "regulation of T cell activation"
## [3] "leukocyte migration"
## [4] "leukocyte mediated cytotoxicity"
## [5] "T cell mediated cytotoxicity"
## [6] "cell killing"
## [7] "positive regulation of T cell activation"
## [8] "leukocyte cell-cell adhesion"
## [9] "positive regulation of leukocyte cell-cell adhesion"
## [10] "positive regulation of leukocyte activation"
2.细胞基因集
自从有了msigdbr这个好东西，就不再需要研究gmt了哦。新加上的C8就是各种细胞的marker基因。
那么难道msigdb的基因集只能做GSEA吗？
那还真不是。
Y叔叔写的函数enrichr，只要你拿到基因与通路/分类名称之间的对应关系，自定义的富集分析，你值得拥有。
C8 <- msigdbr(species = "Homo sapiens",
category = "C8") %>%
dplyr::select(gs_name,gene_symbol )
C8$gs_name = C8$gs_name %>%
str_replace("_","///") %>%
str_split("///",simplify = T) %>%
.[,2] %>%
str_to_lower() %>%
str_replace_all("_"," ")
head(C8)
## # A tibble: 6 × 2
## gs_name gene_symbol
## <chr> <chr>
## 1 liver c10 mvecs 1 A2M
## 2 liver c10 mvecs 1 ACKR1
## 3 liver c10 mvecs 1 ACP5
## 4 liver c10 mvecs 1 ACVRL1
## 5 liver c10 mvecs 1 ADAM15
## 6 liver c10 mvecs 1 ADAMTS1
x <- enricher(g, TERM2GENE = C8)
head(x@result$Description,10)
## [1] "adult olfactory neuroepithelium nk cells"
## [2] "embryonic ctx brain effector t cell"
## [3] "developing heart c9 b t cell"
## [4] "liver c1 nk nkt cells 1"
## [5] "fetal adrenal lymphoid cells"
## [6] "fetal intestine lymphoid cells"
## [7] "fetal placenta lymphoid cells"
## [8] "fetal lung lymphoid cells"
## [9] "adult olfactory neuroepithelium cd8 t cells"
## [10] "fetal stomach lymphoid cells"
dotplot(x)
天高任鸟飞，只要能做出TERM2GENE，就可以实现富集分析自由了。
如果因为代码看不懂，而跟不上正文的节奏，可以来找我，系统学习。我的课程都是循环开课。下一期的时间，点进去咨询微信咯 新的一年，你要不要来找我们学习生信
https://mp.weixin.qq.com/s/hMlLjJha27XSXcMf0hRLOQ