Immugent在推文一文帶你了解單細(xì)胞數(shù)據(jù)基因集打分的所有算法中滓走，介紹了11種常用的對(duì)單細(xì)胞數(shù)據(jù)進(jìn)行基因集打分的算法沛鸵，其中提到了一個(gè)非常給力的R包：irGSEA液斜。它內(nèi)置了"AUCell", "UCell", "singscore", "ssgsea"四種算法，并且可以將多種基因集的富集分?jǐn)?shù)矩陣直接保存到Seurat對(duì)象中，那么這期推文Immugent就給大家用示例數(shù)據(jù)來演示一下這個(gè)包的常見用法蔼两。
首先就是包的安裝問題，可以參考我前一篇文章的代碼。此外算利，Immugent已經(jīng)咨詢過包的作者，最好是用最新版(4.1) 的R上進(jìn)行安裝泳姐，但是小編用R version 4.0.5 (2021-03-31)也裝上這個(gè)包了效拭，只是在調(diào)用其中一個(gè)函數(shù)時(shí)出現(xiàn)了以下問題：
原因是這個(gè)包鏈接到的"GSVA"包函數(shù)是最新版的1.40.1；而4.05版本的R只能裝低版本的“GSVA"包，能裝上的最高版本是1.38.2的缎患。
所以在低版本上即使能裝上irGSEA包慕的，其中的 “ssgsea” 算法是不能使用的，但是從專業(yè)角度看挤渔，"AUCell", "UCell", "singscore"更適用于單細(xì)胞數(shù)據(jù)肮街，所以使用前3種算法就足夠使用的了。
好了判导，下面開始展示具體的用法嫉父，為了每位伙伴都能很好的復(fù)現(xiàn)，小編就用Seuratdata的數(shù)據(jù)進(jìn)行演示眼刃。
-------------------------------1.準(zhǔn)備數(shù)據(jù)---------------------------
library(SeuratData)# view all available datasets
View(AvailableData())
InstallData("pbmc3k")# download 3k PBMCs from 10X Genomics
library(Seurat)# loading dataset
data("pbmc3k.final")
pbmc3k.final <- UpdateSeuratObject(pbmc3k.final)
pbmc3k.final
DimPlot(pbmc3k.final, reduction = "umap",
group.by = "seurat_annotations",label = T) + NoLegend()
我們可以看到這是來源于人外周血的單細(xì)胞數(shù)據(jù)绕辖，包含各種免疫細(xì)胞，其主要分為三大塊：B細(xì)胞鸟整，T細(xì)胞和髓系淋巴細(xì)胞引镊。下面就是要使用上述的四種方法計(jì)算基因集的富集分?jǐn)?shù)了，這里我們使用人的Hallmark數(shù)據(jù)集來進(jìn)行展示篮条。
-------------------------------2.計(jì)算富集分?jǐn)?shù)---------------------------
library(UCell)
library(irGSEA)
Idents(pbmc3k.final) <- pbmc3k.final$seurat_annotations #定義要分析的細(xì)胞亞群
pbmc3k.final <- irGSEA.score(object = pbmc3k.final, assay = "RNA",
slot = "data", seeds = 123, ncores = 1,
min.cells = 3, min.feature = 0,
custom = F, geneset = NULL, msigdb = T,
species = "Homo sapiens", category = "H",
subcategory = NULL, geneid = "symbol",
method = c("AUCell", "UCell", "singscore",
"ssgsea"),
aucell.MaxRank = NULL, ucell.MaxRank = NULL,
kcdf = 'Gaussian')
Seurat::Assays(pbmc3k.final)
這個(gè)包最方便的就是可以通過一個(gè)函數(shù)直接使用4種算法弟头，同時(shí)計(jì)算目的基因集的富集分?jǐn)?shù)，而且函數(shù)還被作者進(jìn)行優(yōu)化涉茧，使得大大節(jié)約了所需計(jì)算時(shí)間赴恨。此外，因?yàn)閮?nèi)置了”msigdb“包伴栓，我們可以方便的計(jì)算人和其它模式物種的各種常見基因集伦连，如GO, KEGG, Hallmark等。
在計(jì)算出所有細(xì)胞的通路富集分?jǐn)?shù)后钳垮，接下來就需要找出差異的通路進(jìn)行重點(diǎn)分析了惑淳。
-------------------------------3.整合差異基因集---------------------------
result.dge <- irGSEA.integrate(object = pbmc3k.final,
group.by = "seurat_annotations",
metadata = NULL, col.name = NULL,
method = c("AUCell","UCell","singscore",
"ssgsea"))
class(result.dge) #> [1] "list"
為了評(píng)估基因集是否在某個(gè)細(xì)胞亞群中富集，irGSEA通過Wilcox檢驗(yàn)計(jì)算富集分?jǐn)?shù)矩陣中的差異基因集（差異基因的過濾標(biāo)準(zhǔn)為為校正后P值小于0.05）饺窿。然后歧焦，又通過RobustRankAggreg包中的秩聚合算法(RRA)對(duì)差異分析的結(jié)果進(jìn)行綜合評(píng)估，并篩選出在大部分基因集富集分析方法中顯著富集的基因集（綜合評(píng)估的過濾標(biāo)準(zhǔn)為P值小于0.05）肚医。但是小編感覺這里只能對(duì)所有的細(xì)胞亞群同時(shí)進(jìn)行統(tǒng)計(jì)學(xué)檢驗(yàn)绢馍，而實(shí)際應(yīng)用中我們往往只關(guān)注某兩群細(xì)胞，因此如果未來能開發(fā)出像Seurat那樣可以做到指定兩群細(xì)胞進(jìn)行統(tǒng)計(jì)學(xué)差異檢驗(yàn)?zāi)蔷透昝懒恕?/p>
這個(gè)包的另一個(gè)亮點(diǎn)就是可以做出很多好看的圖肠套，下面就簡單做出幾個(gè)圖展示一下
-------------------------------4.可視化展示---------------------------
irGSEA.barplot.plot <- irGSEA.barplot(object = result.dge,
method = c("AUCell", "UCell", "singscore",
"ssgsea"))
irGSEA.barplot
halfvlnplot <- irGSEA.halfvlnplot(object = pbmc3k.final,
method = "UCell",
show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
halfvlnplot
ridgeplot <- irGSEA.ridgeplot(object = pbmc3k.final,
method = "UCell",
show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")
ridgeplot
上面說到我們可以對(duì)經(jīng)典數(shù)據(jù)集進(jìn)行分析舰涌，但現(xiàn)實(shí)操作中我們往往想要對(duì)自己構(gòu)建的目的基因集來進(jìn)行富集分析，那么下面小編就進(jìn)行展示如果構(gòu)建自己的目的基因集你稚；此外瓷耙，值得一提的是"UCell", "singscore"函數(shù)可以通過正反向基因同時(shí)進(jìn)行計(jì)算富集分?jǐn)?shù)朱躺，那我們接下來就來看一下效果如何
-------------------------------5.個(gè)性化分析---------------------------
markers <- list()
markers$TcellN <- c("CD3E+","CD3G+","ITGAE-", "ITGAM-")
markers$TcellY <- c("CD3E+","CD3G+","CD19-","CD79A-","CD79B-","FGFBP2-", "KLRF1-","ITGAE-", "ITGAM-")
pbmc3k.final <- irGSEA.score(object = pbmc3k.final, assay = "RNA", slot = "data", seeds = 123, ncores = 1,msigdb=F, custom = T, geneset = markers, method = c( "UCell", "singscore"), kcdf = 'Gaussian')
scatterplot <- irGSEA.density.scatterplot(object = pbmc3k.final,
method = "singscore",
show.geneset = c("TcellN","TcellY")
reduction = "umap")
scatterplot
咋一看用這兩種不同的基因集做出的富集分?jǐn)?shù)很類似，但是細(xì)細(xì)看來就有不一樣的發(fā)現(xiàn)哺徊。如右邊那群髓系淋巴細(xì)胞室琢，在沒有排除B細(xì)胞的marker基因時(shí)的得分低于右側(cè)排除后的，這是因?yàn)樗柘盗馨图?xì)胞不表達(dá)B細(xì)胞的marker落追，用這種方式做出的綜合評(píng)分會(huì)適當(dāng)加分盈滴，而且這種影響隨著基因數(shù)目增加而增加；此外轿钠，雖然T細(xì)胞群和其它群比較顏色上沒有差異巢钓，但是打分范圍右側(cè)高于左側(cè)，也可以認(rèn)為疗垛，髓系的打分沒有變症汹，是B細(xì)胞群變得更低了，這其實(shí)是更有利于發(fā)現(xiàn)差異的贷腕。
好了背镇，說到這，關(guān)于單細(xì)胞基因集打分的故事就說完了泽裳。接下來小編將會(huì)介紹另一個(gè)基于bulk數(shù)據(jù)進(jìn)行打分的強(qiáng)大R包：IOBR瞒斩，敬請(qǐng)期待！
