Immugent在推文一文帶你了解單細(xì)胞數(shù)據(jù)基因集打分的所有算法中，介紹了11種常用的對單細(xì)胞數(shù)據(jù)進(jìn)行基因集打分的算法，其中提到了一個非常給力的R包：irGSEA。它內(nèi)置了"AUCell", "UCell", "singscore", "ssgsea"四種算法，并且可以將多種基因集的富集分?jǐn)?shù)矩陣直接保存到Seurat對象中，那么這期推文Immugent就給大家用示例數(shù)據(jù)來演示一下這個包的常見用法。

首先就是包的安裝問題，可以參考我前一篇文章的代碼。此外，Immugent已經(jīng)咨詢過包的作者，最好是用最新版(4.1) 的R上進(jìn)行安裝，但是小編用R version 4.0.5 (2021-03-31)也裝上這個包了，只是在調(diào)用其中一個函數(shù)時出現(xiàn)了以下問題：

圖片

原因是這個包鏈接到的"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ǐn)?shù)了，這里我們使用人的Hallmark數(shù)據(jù)集來進(jìn)行展示。

-------------------------------2.計算富集分?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)

這個包最方便的就是可以通過一個函數(shù)直接使用4種算法，同時計算目的基因集的富集分?jǐn)?shù)，而且函數(shù)還被作者進(jìn)行優(yōu)化，使得大大節(jié)約了所需計算時間。此外，因為內(nèi)置了”msigdb“包，我們可以方便的計算人和其它模式物種的各種常見基因集，如GO, KEGG, Hallmark等。

在計算出所有細(xì)胞的通路富集分?jǐn)?shù)后，接下來就需要找出差異的通路進(jì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"

為了評估基因集是否在某個細(xì)胞亞群中富集，irGSEA通過Wilcox檢驗計算富集分?jǐn)?shù)矩陣中的差異基因集（差異基因的過濾標(biāo)準(zhǔn)為為校正后P值小于0.05）。然后，又通過RobustRankAggreg包中的秩聚合算法(RRA)對差異分析的結(jié)果進(jìn)行綜合評估，并篩選出在大部分基因集富集分析方法中顯著富集的基因集（綜合評估的過濾標(biāo)準(zhǔn)為P值小于0.05）。但是小編感覺這里只能對所有的細(xì)胞亞群同時進(jìn)行統(tǒng)計學(xué)檢驗，而實際應(yīng)用中我們往往只關(guān)注某兩群細(xì)胞，因此如果未來能開發(fā)出像Seurat那樣可以做到指定兩群細(xì)胞進(jìn)行統(tǒng)計學(xué)差異檢驗?zāi)蔷透昝懒恕?/p>

這個包的另一個亮點就是可以做出很多好看的圖，下面就簡單做出幾個圖展示一下

-------------------------------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

圖片

上面說到我們可以對經(jīng)典數(shù)據(jù)集進(jìn)行分析，但現(xiàn)實操作中我們往往想要對自己構(gòu)建的目的基因集來進(jìn)行富集分析，那么下面小編就進(jìn)行展示如果構(gòu)建自己的目的基因集；此外，值得一提的是"UCell", "singscore"函數(shù)可以通過正反向基因同時進(jìn)行計算富集分?jǐn)?shù)，那我們接下來就來看一下效果如何

-------------------------------5.個性化分析---------------------------

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基因時的得分低于右側(cè)排除后的，這是因為髓系淋巴細(xì)胞不表達(dá)B細(xì)胞的marker，用這種方式做出的綜合評分會適當(dāng)加分，而且這種影響隨著基因數(shù)目增加而增加；此外，雖然T細(xì)胞群和其它群比較顏色上沒有差異，但是打分范圍右側(cè)高于左側(cè)，也可以認(rèn)為，髓系的打分沒有變，是B細(xì)胞群變得更低了，這其實是更有利于發(fā)現(xiàn)差異的。

好了，說到這，關(guān)于單細(xì)胞基因集打分的故事就說完了。接下來小編將會介紹另一個基于bulk數(shù)據(jù)進(jìn)行打分的強(qiáng)大R包：IOBR，敬請期待！