一、單細(xì)胞基因富集分析算法

一個(gè)不算正式的引言：目前來(lái)說(shuō)，基于基因集進(jìn)行分析已經(jīng)開(kāi)發(fā)出來(lái)了很多成體系的R包或者流程，理解來(lái)看仅颇，基因集評(píng)分其實(shí)就是自定義一個(gè)評(píng)分，然后來(lái)衡量目標(biāo)基因集在某組織的表達(dá)情況碘举，進(jìn)而來(lái)推斷其功能富集情況忘瓦，所以說(shuō)，這個(gè)給了我們以提示引颈，算法是一定的耕皮，但是參考基因集可以是不同的，比如說(shuō)鐵死亡蝙场、銅死亡凌停、細(xì)胞衰老等等

目前來(lái)看，常見(jiàn)的針對(duì)單細(xì)胞的基因集富集分析算法有：GSEA售滤、GSVA罚拟、AddModuleScore、AUCell完箩、ssGSEA以及UCell等等赐俗，關(guān)于其中的鑒別以及或多或少的區(qū)別，如果只是單純的跑流程可能并不需要我們掌握弊知。

但是對(duì)于我們分析人員而言阻逮，進(jìn)一步的學(xué)習(xí)可以幫助我們充分的理解算法本身，也可以在有眾多選擇中選擇自己認(rèn)為正確的秩彤，那么接下來(lái)簡(jiǎn)單的來(lái)聊一聊上述這些算法：

1.GSEA

GSEA的輸入文件是一個(gè)表達(dá)矩陣叔扼，我們需要提前把樣本進(jìn)行分組，然后對(duì)所有基因基于分組來(lái)進(jìn)行排序漫雷，那么這個(gè)時(shí)候我們最常使用的排序方式就是通過(guò)LogFC的大小來(lái)進(jìn)行排序（從大到泄细弧），這種排序天生就具有了一種屬性即可以表示兩組間基因的表達(dá)量變化趨勢(shì)珊拼，排序之后的基因列表其頂部可以看作是上調(diào)的差異基因食呻，其底部可以看作是下調(diào)的差異基因

那么這個(gè)時(shí)候GSEA分析其實(shí)回答我們的就是： 我們所選擇的參考基因集是富集在這個(gè)基因排序列表中的頂部還是底部，如果在頂部澎现，那么很顯然我們選擇的參考基因集就是呈現(xiàn)上調(diào)趨勢(shì)仅胞，反之則是下調(diào)趨勢(shì)，一般來(lái)說(shuō)我們選擇的分組都是疾病組和對(duì)照組剑辫，那么有了趨勢(shì)我們就可以回答疾病組表達(dá)的是參考基因集功能（炎癥干旧、衰老、銅死亡）是高表達(dá)還是低表達(dá)妹蔽，這樣其實(shí)就是很好的規(guī)避掉了單純的ORA算法所帶來(lái)的局限性

一椎眯、單細(xì)胞基因富集分析算法

一個(gè)不算正式的引言：目前來(lái)說(shuō)挠将，基于基因集進(jìn)行分析已經(jīng)開(kāi)發(fā)出來(lái)了很多成體系的R包或者流程，理解來(lái)看编整，基因集評(píng)分其實(shí)就是自定義一個(gè)評(píng)分舔稀，然后來(lái)衡量目標(biāo)基因集在某組織的表達(dá)情況，進(jìn)而來(lái)推斷其功能富集情況掌测，所以說(shuō)内贮，這個(gè)給了我們以提示，算法是一定的汞斧，但是參考基因集可以是不同的夜郁，比如說(shuō)鐵死亡、銅死亡粘勒、細(xì)胞衰老等等

目前來(lái)看竞端，常見(jiàn)的針對(duì)單細(xì)胞的基因集富集分析算法有：GSEA、GSVA庙睡、AddModuleScore事富、AUCell、ssGSEA以及UCell等等埃撵，關(guān)于其中的鑒別以及或多或少的區(qū)別赵颅，如果只是單純的跑流程可能并不需要我們掌握虽另。

但是對(duì)于我們分析人員而言暂刘，進(jìn)一步的學(xué)習(xí)可以幫助我們充分的理解算法本身，也可以在有眾多選擇中選擇自己認(rèn)為正確的捂刺，那么接下來(lái)簡(jiǎn)單的來(lái)聊一聊上述這些算法：

1.GSEA

GSEA的輸入文件是一個(gè)表達(dá)矩陣谣拣，我們需要提前把樣本進(jìn)行分組，然后對(duì)所有基因基于分組來(lái)進(jìn)行排序族展，那么這個(gè)時(shí)候我們最常使用的排序方式就是通過(guò)LogFC的大小來(lái)進(jìn)行排序（從大到猩），這種排序天生就具有了一種屬性即可以表示兩組間基因的表達(dá)量變化趨勢(shì)仪缸，排序之后的基因列表其頂部可以看作是上調(diào)的差異基因贵涵，其底部可以看作是下調(diào)的差異基因

那么這個(gè)時(shí)候GSEA分析其實(shí)回答我們的就是： 我們所選擇的參考基因集是富集在這個(gè)基因排序列表中的頂部還是底部，如果在頂部恰画，那么很顯然我們選擇的參考基因集就是呈現(xiàn)上調(diào)趨勢(shì)宾茂，反之則是下調(diào)趨勢(shì)，一般來(lái)說(shuō)我們選擇的分組都是疾病組和對(duì)照組拴还，那么有了趨勢(shì)我們就可以回答疾病組表達(dá)的是參考基因集功能（炎癥跨晴、衰老、銅死亡）是高表達(dá)還是低表達(dá)片林，這樣其實(shí)就是很好的規(guī)避掉了單純的ORA算法所帶來(lái)的局限性

注意：

1.GSEA要求的輸入數(shù)據(jù)是全部表達(dá)矩陣端盆，并不是經(jīng)過(guò)篩選后的表達(dá)矩陣

2.GSEA是一種算法怀骤，如果你是通過(guò)GO數(shù)據(jù)庫(kù)進(jìn)行GSEA富集分析，那么應(yīng)該被叫做GO_GSEA富集分析焕妙，其實(shí)簡(jiǎn)單理解就是蒋伦，我們提供參考基因集然后進(jìn)行GSEA富集分析去探索我們的參考基因集在實(shí)驗(yàn)組究竟是上調(diào)還是下調(diào)

3.對(duì)于是否進(jìn)行批次效應(yīng)的去除，因?yàn)槲覀冃枰M(jìn)行差異分析獲得可以排序的依據(jù)焚鹊，所以這里筆者認(rèn)為我們需要進(jìn)行批次效應(yīng)的去除凉敲，因?yàn)樾枰M(jìn)行差異分析獲得LogFC，但是簡(jiǎn)單來(lái)說(shuō)就是輸入數(shù)據(jù)其實(shí)就是可以進(jìn)行差異分析的輸入數(shù)據(jù)即可

2.GSVA

GSVA富集分析全稱(chēng)叫做基因集變異分析寺旺，在這里說(shuō)一下個(gè)人的簡(jiǎn)單理解爷抓，其實(shí)我們可以把其看作是一個(gè)特殊的富集分析，為什么說(shuō)是特殊呢阻塑，GSVA的觀測(cè)不再是Gene symbol而是轉(zhuǎn)換成了通路名稱(chēng)蓝撇，所以我們可以通過(guò)GSVA獲得表達(dá)矩陣轉(zhuǎn)換成通路表達(dá)矩陣后的結(jié)果，然后我們對(duì)轉(zhuǎn)換后的結(jié)果進(jìn)行差異分析陈莽，就可以得到通路LogFC的情況

注意：

1.GSVA選擇不同的參考基因集得到的結(jié)果是不一樣的渤昌，這一點(diǎn)我們可以靈活運(yùn)用

2.輸入數(shù)據(jù)可以是counts也可以是logTPM，如果是counts則參數(shù)kcdf需要選擇"Poisson"走搁，如果是logTPM則參數(shù)kcdf需要選擇"Gaussian"

3.AddModuleScore

關(guān)于AddModuleScore函數(shù)独柑，我曾經(jīng)寫(xiě)過(guò)系統(tǒng)的推文，感興趣的小伙伴可以直接點(diǎn)擊下方鏈接進(jìn)行跳轉(zhuǎn)閱讀：

Seurat包的打分函數(shù)AddModuleScore

雖然可能理解算法流程對(duì)我們使用可能沒(méi)有太多的幫助私植，但是這里依舊想要簡(jiǎn)單的概括一下算法的大概流程：通過(guò)計(jì)算參考基因集在表達(dá)矩陣中所有基因的平均表達(dá)量忌栅，然后通過(guò)平均表達(dá)量把表達(dá)矩陣切割成若干份，然后從切割的部分中隨機(jī)抽取這些切塊作為背景基因曲稼，知道這些的重點(diǎn)其實(shí)是下面這句話(huà)既然其中存在隨機(jī)索绪，那么即使使用相同參考基因集為相同數(shù)據(jù)進(jìn)行打分，也可能會(huì)產(chǎn)生不同的結(jié)果

4.AUCell

首先我們通過(guò)盡量簡(jiǎn)單的語(yǔ)言來(lái)回答一下AUCell算法贫悄，其實(shí)本質(zhì)上就是計(jì)算我們所選擇的參考基因集是否在單細(xì)胞中每個(gè)變量上富集瑞驱，因?yàn)閱渭?xì)胞數(shù)據(jù)的變量是細(xì)胞ID，我們也可以說(shuō)是計(jì)算是否在細(xì)胞上富集

簡(jiǎn)單來(lái)說(shuō)AUCell算法分為三個(gè)步驟：

1.針對(duì)每個(gè)細(xì)胞中的每個(gè)基因進(jìn)行排序（既然涉及到了排序窄坦，那么顯然標(biāo)準(zhǔn)化與否對(duì)于我們的影響就可以忽略不計(jì)了唤反，因?yàn)榧词箻?biāo)準(zhǔn)化也不會(huì)改變誰(shuí)大誰(shuí)小這個(gè)直觀的事實(shí)）

2.用過(guò)使用AUC來(lái)評(píng)估我們的參考基因集是否在每個(gè)細(xì)胞前5%表達(dá)基因集內(nèi)富集，然后進(jìn)而來(lái)評(píng)估參考基因集針對(duì)每個(gè)細(xì)胞的富集情況

5.ssGSEA

ssGSEA全名其實(shí)被稱(chēng)為"單樣本基因集富集分析",從名字其實(shí)我們就可以看出鸭津，其是GSEA方法的擴(kuò)展彤侍，既然是擴(kuò)展，那么我們就首先需要知道GSEA算法的局限性

GSEA算法的局限其實(shí)就是我們需要提供排序信息曙博，而我們的排序信息在很多情況下其實(shí)源自不同分組下差異表達(dá)拥刻，也就是logFC，那么這個(gè)時(shí)候問(wèn)題來(lái)了：如果我只有一個(gè)分組父泳，那么我如何評(píng)估某個(gè)基因集在單樣本（一個(gè)分組）中的表達(dá)情況呢般哼？

那么這個(gè)時(shí)候自然而然就誕生了ssGSEA吴汪，該算法可以讓單個(gè)樣本獲得對(duì)應(yīng)參考基因集的富集評(píng)分，通過(guò)這種形式蒸眠，我們就可以把單一樣本的基因表達(dá)譜轉(zhuǎn)換為基因富集分?jǐn)?shù)表達(dá)譜漾橙，那么這個(gè)時(shí)候我們就可以直觀描述單一樣本下針對(duì)參考基因集富集分?jǐn)?shù)的高低，那么這個(gè)時(shí)候問(wèn)題又來(lái)了：如果我有兩個(gè)分組楞卡，我使用ssGSEA霜运，可以嗎？

那么這個(gè)時(shí)候我們可以參考免疫浸潤(rùn)分析答案自然就呼之欲出了蒋腮，也就是說(shuō)我如果把參考基因集換成免疫相關(guān)基因淘捡，然后通過(guò)使用ssGSEA針對(duì)疾病組和對(duì)照組，那么我們就可以得到疾病組和對(duì)照組分別的免疫浸潤(rùn)情況池摧，然后通過(guò)進(jìn)行秩和檢驗(yàn)焦除，自然就可以得到下面這張?jiān)谖墨I(xiàn)中經(jīng)常出現(xiàn)的小提琴圖：

6.UCell

看到這里，相信我們都或多或少會(huì)有一種感覺(jué)作彤，算法很多膘魄，但是都是得到一個(gè)打分，一個(gè)可以評(píng)估我們參考基因集在樣本中表達(dá)情況的打分竭讳，那么很顯然有了這層理解后创葡，我們后面哪怕再接觸其它"打分"算法，我們也可以知道我們的重心更應(yīng)該放在參考基因集的選擇上绢慢，也就是我們想要明確哪些功能在我們的樣本中

那么回到算法本身灿渴，UCell算法具有以下特點(diǎn)：

1.對(duì)數(shù)據(jù)集本身要求不高，并不會(huì)受到數(shù)據(jù)集異質(zhì)性的影響

2.對(duì)于處理大型數(shù)據(jù)具有一定優(yōu)勢(shì)

3.可以方便銜接Seurat包呐芥，快捷方便

好逻杖，那么介紹到這里，相信各位小伙伴依然還有疑惑思瘟，上述看似干貨的東西顯然并不能使我們掌握單細(xì)胞基因集富集分析，那么這里晨曦自然并不是僅僅止步于此

二闻伶、R包實(shí)踐

接下來(lái)給大家介紹一個(gè)可以自用定制單細(xì)胞基因集富集分析需求同時(shí)可視化也很美觀的R包——irGSEA包

接下來(lái)將直接演示代碼滨攻，

library(Seurat)
library(SeuratData)
library(UCell)
library(irGSEA)
# loading datasetdata( "pbmc3k.final") 
pbmc3k. final<- UpdateSeuratObject(pbmc3k. final) # 
plotDimPlot(pbmc3k. final, reduction = "umap", group.by = "seurat_annotations",label = T) + NoLegend 

這個(gè)數(shù)據(jù)就是Seurat官網(wǎng)上經(jīng)典的PBMC數(shù)據(jù)，我們這里演示的流程可以無(wú)縫對(duì)接自己的Seurat流程數(shù)據(jù)

Idents(pbmc3k. final) <- pbmc3k. final$seurat_annotations 
#設(shè)置亞群蓝翰，后續(xù)選擇注釋后的數(shù)據(jù)
#計(jì)算評(píng)分
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') 
#整合差異基因集
result.dge <- irGSEA.integrate(object = pbmc3k. final, group.by = "seurat_annotations", metadata = NULL, col.name = NULL, method = c( "AUCell", "UCell", "singscore", "ssgsea"))

一步函數(shù)即可完成相關(guān)分析光绕，在這里作者內(nèi)置了MSigDB數(shù)據(jù)庫(kù)的相關(guān)基因集，當(dāng)然R包本身也支持自定義基因集畜份，這一塊就給了我們無(wú)限的可能诞帐，其實(shí)就是更加的自由

#可視化
irGSEA.heatmap.plot < -irGSEA.heatmap( object= result.dge,method= "RRA", top= 50,show.geneset= NULL)
irGSEA.heatmap.plot

#可視化
irGSEA.barplot.plot <- irGSEA.barplot( object= result.dge, method = c( "AUCell", "UCell", "singscore", "ssgsea")) 
irGSEA.barplot.plot

當(dāng)然還有更多可視化的展示形式，大家感興趣可以去作者的github上進(jìn)行進(jìn)一步的學(xué)習(xí)......

生活很好爆雹，有你更好