R語言實現(xiàn)通路富集打分

terminator_523 2019-12-02

展開全文

我們大家應(yīng)該對通路富集分析都很熟悉，比如GSEA,，DAVID等,。都是在大量文章中常見的通路富集方法，那么今天我們也給大家介紹一個更加復(fù)雜的通路富集分析的前期數(shù)據(jù)處理包GSVA（gene set variation analysis）,。是一種非參數(shù)的無監(jiān)督分析方法,，主要用來評估芯片核轉(zhuǎn)錄組的基因集富集結(jié)果。主要是通過將基因在不同樣品間的表達量矩陣轉(zhuǎn)化成基因集在樣品間的表達量矩陣,，從而來評估不同的通路在不同樣品間是否富集,。具體的一個分析流程如下：

首先我們看下安裝，在R語言3.5版本以上的安裝代碼如下：

if(!requireNamespace('BiocManager', quietly = TRUE)) install.packages('BiocManager') BiocManager::install('GSVA')BiocManager::install('GSVAdata')

接下來我們看下前期的數(shù)據(jù)清洗,。我們需要用到R包genefilter中的nsFilter函數(shù)：

其中主要的參數(shù)：

Require.GOBP, require.GOCC, require.GOMF, require.CytoBand指是否只保留具有相對應(yīng)功能的基因,，默認是FALSE。

remove.dupEntrez 移除多個探針指向一個基因的情況,。移除的標準是var.func后面的計算方法包括：overall mean, median,variance, IQR,。

Var.cutoff 指需要排除的不穩(wěn)定的探針部分。

Var.filter 邏輯值,，如果是TRUE就會調(diào)用var.func進行計算,。

FilterByQuantile 邏輯值確定對var.cutoff進行篩選計算。

接下來我們看下數(shù)據(jù)的前期預(yù)處理：

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')##此處默認值是模型認為最好的一種設(shè)置,。

leukemia_filtered_eset <-filtered_eset$eset

接下來我們看下我們要用到的核心函數(shù)gsva：

其中主要的參數(shù)：

Gset.idx.list 自己需要的基因集列表,，這個可以自己進行定義，主要是GSEA提供的數(shù)據(jù)對象，主要數(shù)據(jù)來源GSEAbase包,。

首先載入數(shù)據(jù)：

data(c2BroadSets)

然后是對通路數(shù)據(jù)進行篩選：

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網(wǎng)站進行獲?。?/p>

geneSets <-getGmt('test.geneset')

Method 在實現(xiàn)GSVA的同時還實現(xiàn)了其他相關(guān)的計算方法，我們就不一一介紹了,。

Kcdf 指的是數(shù)據(jù)類型選擇,。RNA-seq的原始整數(shù)的read count 在使用gsva時需要設(shè)置kcdf='Possion',如果是取過log的RPKM,TPM等結(jié)果可以使用默認的值。

接下來就是算法的計算過程,，實例如下：

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)