以下文章來源于生信會客廳 ，作者Kinesin
生信會客廳分享單細(xì)胞轉(zhuǎn)錄組、空間轉(zhuǎn)錄組、轉(zhuǎn)錄調(diào)控及微生物多樣性的分析流程、代碼及相關(guān)生信知識。
上期專題我們介紹了單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的基礎(chǔ)分析，然而那些分析只是揭開了組織異質(zhì)性的面紗，還有更多的生命奧秘隱藏在數(shù)據(jù)中等待我們發(fā)掘。本專題將介紹一些單細(xì)胞轉(zhuǎn)錄組的高級分析內(nèi)容：多樣本批次校正、轉(zhuǎn)錄因子分析、細(xì)胞通訊分析、基因集變異分析和更全面的基因集富集分析。不足之處請大家批評指正，歡迎添加Kinesin微信交流探討!
GSEA與GSVA簡介
**基因集的概念 **
GSEA全稱Gene Set Enrichment Analysis，GSVA全稱Gene Set Variation Analysis，它們都是基于基因集開展的分析，因此我們先要了解基因集的定義?；蚣櫭剂x就是一些基因的集合，任何一些基因放在一起都可以叫做基因集，但是我們用來分析的基因集要求有一定的生物學(xué)意義。最常見的基因集數(shù)據(jù)庫如GO與KEGG，它們一個按照基因本體論將基因分門別類，一個按照代謝通路將相關(guān)基因集合在一起。除此之外我們還可以按轉(zhuǎn)錄因子調(diào)控網(wǎng)絡(luò)、共表達(dá)網(wǎng)絡(luò)、定義生物狀態(tài)的marker基因列表等把基因集合成有一定生物學(xué)意義的基因集。
**MSigDB基因集數(shù)據(jù)庫 **
GSEA是由Broad研究所開發(fā)的一種富集方法，他們在提出該方法的同時還提供了一個基因集數(shù)據(jù)庫——MSigdb。它從位置，功能，代謝途徑，靶標(biāo)結(jié)合等多種角度出發(fā)，構(gòu)建出了許多的基因集合，Broad研究所將他們構(gòu)建的基因集合保存在MSigDB，官網(wǎng)地址如下：http://software.broadinstitute.org/gsea/msigdb/index.jsp在MSigDB中，將所有的基因集劃分為以下8大類別：
1. H：hallmark gene sets
特征基因集，由定義生物狀態(tài)和進(jìn)程的marker基因組成。
2. C1：positional gene sets
位置基因集，包含人類每條染色體上的不同cytoband區(qū)域?qū)?yīng)的基因集合。
3. C2：curated gene sets
代謝通路基因集，包含KEGG, Reactome, BioCarta數(shù)據(jù)庫，以及文獻(xiàn)和專家支持的基因集信息。
4. C3：motif gene sets
靶基因集，包含了miRNA靶基因集和轉(zhuǎn)錄因子調(diào)控基因集兩大類。
5. C4：computational gene sets
計算基因集，計算機軟件預(yù)測出來的基因集，主要是和癌癥相關(guān)的基因。
6. C5：GO gene sets
基因本體基因集，包含了Gene Ontology對應(yīng)的基因集合。
7. C6：oncogenic signatures
癌癥擾動基因集，來源于藥物處理腫瘤后基因差異表達(dá)數(shù)據(jù)，包含已知條件處理后基因表達(dá)量發(fā)生變化的基因。
8. C7：immunologic signatures
免疫基因集，包含了免疫系統(tǒng)功能相關(guān)的基因集合。
**GSEA的分析原理 **
常規(guī)GO/KEGG富集分析需要設(shè)定閾值過濾差異基因，閾值太寬富集的結(jié)果太多，閾值太嚴(yán)又可能會遺漏一些關(guān)鍵結(jié)果。GO/KEGG富集的結(jié)果通常還很寬泛，并不能很好地解釋生物學(xué)現(xiàn)象。有鑒于此，Broad研究所開發(fā)了基因集富集分析(GSEA)方法。GSEA使用無監(jiān)督算法，不用過濾任何基因，配合MSigDB數(shù)據(jù)庫使用，更容易找到解釋生物學(xué)現(xiàn)象的基因集。其原理如下：
GSEA分析要先將樣本做組間對比分析，GSEA自帶9種分析方法，分為基因表達(dá)值差異分析和相關(guān)性分析兩大類。對于Case/Contral的實驗設(shè)計，差異分析方法更為常用，這其中又以默認(rèn)的信噪比和大家熟悉的差異倍數(shù)用的最多。對比分析之后要按結(jié)果將基因排序，以差異倍數(shù)方法為例，把所有基因按差異倍數(shù)(FC)的值降序排列以供后續(xù)分析。上圖小人腳下的小方塊代表排序好的差異基因列表，藍(lán)色之外的其他色塊代表屬于某個基因集的基因，如黃色屬于基因集A，綠色屬于基因集B。最下面高低不等的豎條代表與基因列表對應(yīng)的FC值，紅色上調(diào)、藍(lán)色下調(diào)、黃色沒有變化?；蚣母患治鲂枰?jīng)歷三步：
基因集A富集分析時，小人從基因列表的左端走到右端，每經(jīng)過一個藍(lán)色基因扣分，每遇到一個黃色基因加分，扣分時與FC無關(guān)，加分時考慮FC的權(quán)重?；蚣疉最終的富集分?jǐn)?shù)(ES)是小人曾經(jīng)得過的最高/低分，實際公式比這復(fù)雜，但基本理念如此。
采用置換檢驗計算基因集A的顯著性，即p值。
基因集A富集分析完成后，按上述同樣的方法完成基因集B、C直至所有輸入基因集的分析。所有需要富集分析的基因集都計算ES和p值之后，將ES轉(zhuǎn)換為標(biāo)準(zhǔn)富集分?jǐn)?shù)(NES)，并計算校正后p值。
聽完我的解釋之后再看官方的解釋可能更容易理解：
A GSEA overview illustrating the method. (A) An expression dataset sorted by correlation with phenotype, the corresponding heat map, and the ‘‘gene tags,’’ i.e., location of genes from a set S within the sorted list. (B) Plot of the running sum for S in the dataset, including the location of the maximum enrichment score (ES) and the leading-edge subset.
**GSVA的分析原理 **
GSEA雖然是一種強大的富集分析工具，但是它的應(yīng)用場景通常局限于Case/Control的實驗設(shè)計。對于表型（分組）復(fù)雜的大樣本量研究，例如TCGA和單細(xì)胞圖譜這樣的項目，分析起來就非常困難。因此，Broad研究所在GSEA發(fā)布8年之后，又開發(fā)了GSVA算法來拓展基因集分析的應(yīng)用。GSVA不需要預(yù)先進(jìn)行樣本之間的差異分析，它依據(jù)表達(dá)矩陣就可以計算每個樣本中特定基因集的變異分?jǐn)?shù)。簡單的說，輸入以基因為行的表達(dá)矩陣和基因集數(shù)據(jù)庫給GSVA，它就輸出以基因集名稱為行的變異分?jǐn)?shù)矩陣，如下圖所示：
左側(cè)輸入基因表達(dá)矩陣和基因集數(shù)據(jù)庫，中間是GSVA算法原理，右側(cè)是輸出的基因集變異分?jǐn)?shù)矩陣?；蚣儺惙?jǐn)?shù)可以理解為基因集內(nèi)所有基因的綜合表達(dá)值。
GSEA的安裝與使用
GSEA官方提供很多版本，但是不提倡R語言版，下載網(wǎng)址如下：
https://www.gsea-msigdb.org/gsea/downloads.jsp
**準(zhǔn)備分析文件 **單細(xì)胞GSEA分析的難點是輸入文件的準(zhǔn)備，大家在R語言中運行下列代碼，就可以得到符合GSEA要求的特定格式文件。
##R語言準(zhǔn)備gsea輸入文件
準(zhǔn)備好的文件格式：表達(dá)數(shù)據(jù)gct文件內(nèi)容節(jié)選
分組cls文件內(nèi)容節(jié)選
得到上面的gct和cls文件，再加上基因集數(shù)據(jù)庫gmt文件，就可以進(jìn)行GSEA分析了。建議在java版本上運行，官方的R包太簡陋了，java版上的好多功能都不支持。
** R語言版代碼與結(jié)果展示 **
##安裝官方R包GSEA
分析得到的結(jié)果很多，分為結(jié)果總覽，分組概況，以及各個基因集的富集情況詳細(xì)。這里展示NES分?jǐn)?shù)最高的通路的富集結(jié)果：
上圖對應(yīng)數(shù)據(jù)
**java版分析代碼 **
提前安裝好java11，下載圖中所示GSEA版本，解壓即可
#以下命令在linux環(huán)境中運行，windows環(huán)境用gsea-cli.bat代替gsea-cli.sh
**圖形界面版GSEA操作 **
表達(dá)數(shù)據(jù)和分組文件準(zhǔn)備好之后，圖形界面版操作沒有難度，簡單介紹一下：
這是開始界面，先加載數(shù)據(jù)
點擊箭頭所指按鈕之后，就會彈出對話框讓你選擇文件，分次輸入我們之前準(zhǔn)備好的表達(dá)數(shù)據(jù)文件(.gct)和分組文件(.cls)
選好文件之后，文件目錄會出現(xiàn)在左右兩邊箭頭所指的地方。如果右邊沒有，雙擊左邊即可。我這兒左邊有三個文件是因為加載了本地保存的基因集數(shù)據(jù)庫文件。
文件加載完畢之后，點擊箭頭所指運行會出現(xiàn)右邊的對話框，注意一下紅框標(biāo)注的選項。Collapse/Remap項選擇“No_Collapse”（Collapse適用芯片數(shù)據(jù)），Permutation項選擇“phenotype”（更適合單細(xì)胞數(shù)據(jù)），Chip項可不選或選擇有“Human_Symbol_with_”關(guān)鍵字的項目（選擇之后結(jié)果中基因有簡短注釋）。
Basic fields點開之后選擇保存結(jié)果的目錄，點擊“Run”就開始運行了。java版（圖形界面的也是java版）運行之后結(jié)果與R語言版一致，不過分析報告是網(wǎng)頁版，點擊超鏈接就可以查看細(xì)節(jié)，使用起來非常方便。
圖形也更美觀
GSVA的安裝與使用
** 安裝GSVA **
#從bioconductor安裝
**準(zhǔn)備輸入數(shù)據(jù) **
GSVA要求輸入表達(dá)矩陣和基因集列表。表達(dá)矩陣從seurat對象導(dǎo)入即可，以列表的形式提供基因集數(shù)據(jù)有一定難度，我編寫了一個函數(shù)解決此問題。
library(Seurat)
**結(jié)果可視化 **
GSVA的分析結(jié)果相當(dāng)于表達(dá)矩陣，用熱圖來可視化比較常用，為了體現(xiàn)單細(xì)胞的特點我也用FeaturePlot展示了幾個基因集。
library(pheatmap)
隨意挑了幾個基因集展示gsva富集分?jǐn)?shù)
GSVA分析給我們的不是階段性結(jié)果，而是觀察數(shù)據(jù)的一個新的視角，讓我們以基因集為單位研究數(shù)據(jù)?；诨蚣母患仃?，我們可以做差異分析、表型相關(guān)性分析、生存分析，也可以基于基因集的功能來鑒定細(xì)胞類型。
References
https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideFrame.htmlhttps://www.gsea-msigdb.org/gsea/msigdb/collection_details.jspLaurila E , Mootha V K , Lindgren C M , et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes[J]. 2003, 34(3):267-273.Subramanian, Tamayo, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. 2005, PNAS 102, 15545-15550.Hnzelmann S , Castelo R , Guinney J . GSVA: gene set variation analysis for microarray and RNA-Seq data[J]. Bmc Bioinformatics, 2013, 14(1):7-7.
