文章發(fā)表于2013年,，GSVA: gene set variation analysis for microarray and RNA-Seq data 同樣是broad 研究生出品,，其在2005年PNAS發(fā)表的gsea已經(jīng)高達1.4萬的引用了，不過這個GSVA才不到300,。
算法細節(jié) 算法本身就不是很好理解,，并不強求一定要理解透徹，可以參考2005年的GSEA算法：
GSVA starts by evaluating whether a gene i is highly or lowly expressed in sample j in the context of the sample population distribution.
可以是芯片雜交的信號代表的表達量,，也可以是轉錄組測序定量,。
For each gene expression profile xi={xi1,…,xin}, a non-parametric kernel estimation of its cumulative density function is performed.
We offer two approaches for turning the KS like random walk statistic into an enrichment statistic (ES) (also called GSVA score), the classical maximum deviation method and a normalized ES.
而且作者也在測試數(shù)據(jù)和真實數(shù)據(jù)把自己的GSVA算法跟GSEA，PLAGE, single sample GSEA (ssGSEA) 或者其它算法進行了比較,， 還在TCGA的ovarian serous cystadenocarcinoma (OV)癌癥表達矩陣(n =588) ,，用MSigDB數(shù)據(jù)庫的 canonical gene sets (C2) 基因集做了比較和測試。
還比較了轉錄組測序數(shù)據(jù)和芯片數(shù)據(jù),，這些數(shù)據(jù)都提供了下載鏈接,，最后作者把算法打包成了 Bioconductor package for R under the name GSVA at http://www..
安裝GSVA這個R包 安裝并且查看21頁的PDF教程：
## try http:// if https:// URLs are not supported最新版教程：https://www./packages/devel/bioc/vignettes/GSVA/inst/doc/GSVA.pdf
其實核心函數(shù)就是gsva()，需要兩個輸入：the gene expression data and a collection of gene sets.
其實這個函數(shù)也可以選擇其它3個模型：
method="plage" (Tomfohr et al., 2005). Pathway level analysis of gene expression (PLAGE)
method="zscore" (Lee et al., 2008). The combined z-score method a
method="ssgsea" (Barbie et al., 2009). Single sample GSEA (ssGSEA) calculates a gene set
另外一個比較重要的參數(shù)是： default argument mx.diff=TRUE to obtain approximately
非常多的文章都在引用該算法，比如：https://www./articles/srep16238#f1
先在模擬數(shù)據(jù)應用GSVA 代碼很簡單,，構造一個 30個樣本,，2萬個基因的表達矩陣,， 加上 100 個假定的基因集。
library(GSVA)這樣就可以檢驗我們假定的100個基因集在我們的30個樣本的GSVA score值分布情況,。
值得注意的是,，這里的gsva函數(shù)接受的是一個純粹的表達矩陣matrix和一個純粹基因集合list，實際上通常是一個 ExpressionSet 和 GeneSetCollection 對象,，所以大家務必學會R里面的對象哈,，不然永遠只能看懂簡單代碼。
兩個模型的區(qū)別 par(mfrow=c(1,2), mar=c(4, 4, 4, 1))在真實數(shù)據(jù)(白血病芯片數(shù)據(jù))上面測試GSVA算法 library(GSEABase)是17個MLL病人和20個ALL病人的古老的affymetrix芯片表達數(shù)據(jù),，共12626個探針，這樣的表達矩陣首先需要過濾哈,，這里直接用 genefilter 包提供的 nsFilter函數(shù)來進行過濾,。
filtered_eset <- nsFilter(leukemia_eset, require.entrez=TRUE, 剩下的是 4292個探針的表達矩陣。
根據(jù)表型數(shù)據(jù)使用limma包來找到有顯著差異的基因集 因為每個基因集都在每個樣本里面得到了一個值,，所以這時候相當于有了一個新的表達矩陣,，而且這些樣本的表型數(shù)據(jù)仍然是存在的，所以可以借鑒差異分析的算法了,。
cache(leukemia_es <- gsva(leukemia_filtered_eset, c2BroadSets,Thus, there are 38 MSigDB C2 curated pathways that are differentially activated between MLL and ALL
同樣limma當然是可以找顯著差異表達的基因的 logFCcutoff <- log2(2)Here, 122 genes show up as being differentially expressed with a minimum fold-change of 2 at 0.1% FDR.
可以看到,，兩個代碼唯一的變化就是 leukemia_filtered_eset 和 leukemia_es而已。這樣的差異分析結果同樣也是可以畫火山圖,，熱圖,，代碼就不贅述了，非常簡單,。
先看兩個火山圖的區(qū)別：
然后看兩個熱圖的區(qū)別,；
還可以對包內(nèi)置的TCGA數(shù)據(jù)進行測試 data(gbm_VerhaakEtAl)exprs(gbm_eset)[1:4,1:4] ### 如下 assayData: 11861 features, 173 samples protocolData: none featureData: none experimentData: use 'experimentData(object)' Annotation: 同樣的方式做gsva分析
data(gbm_VerhaakEtAl)不同算法在轉錄組測序數(shù)據(jù)的表現(xiàn) 前面我們說到過gsva函數(shù)還提供了另外3個算法,，這里就不細細講解了,。
也可以只關注3個主流通路數(shù)據(jù)庫：
canonicalC2BroadSets <- c2BroadSets[c(grep("^KEGG", names(c2BroadSets)),還可以增加自己感興趣的基因集
data(genderGenesEntrez)也可以去broad官網(wǎng)里面下載最新版gmt文件，然后讀入到R語言
myC7 <- getGmt("c7.all.v5.1.entrez.gmt", geneIdType=EntrezIdentifier(), collectionType=BroadCollection(category="c7"), sep="\t")這個時候需要詳細了解 GSEABase 包的設計,。
library(GSEABase)
