R语言实现通路富集打分
我们大家应该对通路富集分析都很熟悉，比如GSEA，DAVID等。都是在大量文章中常见的通路富集方法，那么今天我们也给大家介绍一个更加复杂的通路富集分析的前期数据处理包GSVA（gene set variation analysis）。是一种非参数的无监督分析方法，主要用来评估芯片核转录组的基因集富集结果。主要是通过将基因在不同样品间的表达量矩阵转化成基因集在样品间的表达量矩阵，从而来评估不同的通路在不同样品间是否富集。具体的一个分析流程如下：
首先我们看下安装，在R语言3.5版本以上的安装代码如下：
if(!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GSVA")
BiocManager::install("GSVAdata")接下来我们看下前期的数据清洗。我们需要用到R包genefilter中的nsFilter函数：
其中主要的参数：
Require.GOBP, require.GOCC, require.GOMF, require.CytoBand指是否只保留具有相对应功能的基因，默认是FALSE。
remove.dupEntrez 移除多个探针指向一个基因的情况。移除的标准是var.func后面的计算方法包括：overall mean, median,variance, IQR。
Var.cutoff 指需要排除的不稳定的探针部分。
Var.filter 逻辑值，如果是TRUE就会调用var.func进行计算。
FilterByQuantile 逻辑值确定对var.cutoff进行筛选计算。
接下来我们看下数据的前期预处理：
library(GSVAdata)
library(GSVA)
data(leukemia)
leukemia_eset
filtered_eset <- nsFilter(leukemia_eset,require.entrez=TRUE, remove.dupEntrez=TRUE,var.func=IQR, var.filter=TRUE,var.cutoff=0.5, filterByQuantile=TRUE,feature.exclude="^AFFX")##此处默认值是模型认为最好的一种设置。
leukemia_filtered_eset <-filtered_eset$eset
接下来我们看下我们要用到的核心函数gsva：
其中主要的参数：
Gset.idx.list 自己需要的基因集列表，这个可以自己进行定义，主要是GSEA提供的数据对象，主要数据来源GSEAbase包。
首先载入数据：
data(c2BroadSets)
然后是对通路数据进行筛选：
canonicalC2BroadSets <-c2BroadSets[c(grep("^KEGG",names(c2BroadSets)),grep("^REACTOME",names(c2BroadSets)),grep("^BIOCARTA", names(c2BroadSets)))]
data(genderGenesEntrez)
MSY <- GeneSet(msYgenesEntrez,geneIdType=EntrezIdentifier(),collectionType=BroadCollection(category="c2"),setName="MSY")
XiE <- GeneSet(XiEgenesEntrez,geneIdType=EntrezIdentifier(),collectionType=BroadCollection(category="c2"),setName="XiE")
canonicalC2BroadSets <-GeneSetCollection(c(canonicalC2BroadSets, MSY, XiE))
另外一种获取基因集的方式是通过GSEA网站进行获取：
geneSets <-getGmt("test.geneset")Method 在实现GSVA的同时还实现了其他相关的计算方法，我们就不一一介绍了。
Kcdf 指的是数据类型选择。RNA-seq的原始整数的read count 在使用gsva时需要设置kcdf="Possion",如果是取过log的RPKM,TPM等结果可以使用默认的值。
接下来就是算法的计算过程，实例如下：
data(commonPickrellHuang)
canonicalC2BroadSets <-c2BroadSets[c(grep("^KEGG", names(c2BroadSets)))]
#使用GSVA方法进行计算
esmicro <-gsva(huangArrayRMAnoBatchCommon_eset, canonicalC2BroadSets, min.sz=5,max.sz=500,mx.diff=TRUE, verbose=FALSE, parallel.sz=1)
esrnaseq <-gsva(pickrellCountsArgonneCQNcommon_eset, canonicalC2BroadSets, min.sz=5,max.sz=500,kcdf="Poisson", mx.diff=TRUE, verbose=FALSE,parallel.sz=1)
这就是我们的数据结果，行名是样本名称，列名是通路名称。然后我们还可以利用pheatmap将数据进行可视化：
pheatmap::pheatmap(esrnaseq)
至此就得到了我们想要的数据结果，后面的使用计算就是表达矩阵的思路了。