GSAR基因集分析Bioconductor.docx

GSAR基因集分析Bioconductor.docx

《GSAR基因集分析Bioconductor.docx》由会员分享，可在线阅读，更多相关《GSAR基因集分析Bioconductor.docx（24页珍藏版）》请在冰豆网上搜索。

GSAR基因集分析Bioconductor

用R语言做基因集分析——GSAR软件包

GeneSetAnalysisinR--theGSARPackage

YasirRahmatallahandGalinaGlazko

DepartmentofBiomedicalInformatics,UniversityofArkansasforMedicalSciences,LittleRock,AR72205.

yrahmatallah@uams.edu,gvglazko@uams.edu

GSAR版本1.10.0最新修改2017-01-03

GSARversion1.10.0（Lastrevision2017-01-03）

翻译：

任重鲁南方医科大学南方医院,

renzhonglu@,2017-07-13

目录

1简介2

2最小生成树（MST）3

2.1FirstMST3

2.2MST2forcorrelationandPPInetworks4

3统计方法6

3.1Wald-Wolfowitz检验6

3.2Kolmogorov-Smirnov检验6

3.3均值偏差检验Meandeviationtest7

3.4凝聚的F检验AggregatedF-test8

3.5基因集网络相关性分析8

3.5.1方法8

3.5.2零标准偏差存在的问题9

4处理RNA测序计数数据的注释9

5个案研究Casestudies10

5.1p53数据集10

5.1.1简介10

5.1.2数据过滤和标准化10

5.1.3GSA10

5.2急性粒细胞白血病数据集（ALL）13

5.2.1简介13

5.2.2数据过滤和标准化13

5.2.3选定的基因集14

5.3Pickrell数据集15

5.3.1简介15

5.3.2数据过滤和标准化16

5.3.3检验选定的通路17

6会话信息（Sessioninfo）18

7.参考文献（References）19

1.简介（Introduction）

本手册提供了R平台下的GSAR软件包的简要介绍，GSAR包提供了一系列对自包含基因集分析（self-containedgenesetanalysis）的多元统计检验方法。

GSAR用两样本多元非参数统计方法来检验一个原假设和特定的备择假设，比如均值（mean（shift））的差异（位移上的差异），方差（variance（scale））的差异（尺度上的差异）或者相关性结构（correlationstructure）的差异。

该软件包也提供图形可视化工具来展示从表达数据得到的相关性网络，基于最小生成树（minimumspanningtrees）的方法来检验同一基因集在两个条件下的网络相关性结构的改变。

这个可视化工具也可以用在蛋白质-蛋白质相互作用（PPI）网络的分析上来标出最有影响力的蛋白。

这个软件包实现了前文[1–3]中提出的方法，这些方法都在前文[1,2]中经过模拟数据和真实的基因芯片数据集的彻底检验。

图1展示了该包的大纲。

包中的检验函数一次分析一个基因集，包装函数TestGeneSets则可以接受基因集列表并且调用相应的检验函数按次序地对所有的基因集进行检验。

只要将RNA-seq计数数据进行适当地标准化（样本内考虑基因长度的不同（genelength），样本间考虑测序数据量的不同（librarysize）），软件包中的方法也可以应用在处理RNA-seq的计数数据的分析上。

然而对于RNA-seq数据来说，方差是表达均值的一个函数并且对RNA-seq数据使用何种方差检验方法（RKStest，RMDtest，AggrFtest）则高度依赖于标准化的方法。

我们不在这里深入地讨论和描述这个问题。

本手册深入浅出地简要介绍统计方法背后的理论概念，然后提供了多个完整的从基因表达芯片到RNA-seq计数数据的案例分析。

图1：

GSAR的大纲。

统计检验的输入可以是

（1）经过标准化的基因表达芯片或RNA-seq数据的单个基因集的基因表达矩阵，每个样本都有标签来表明样本的所属类别（处理组，正常组）；输入也可以是

（2）所有基因的表达矩阵，每个样本拥有标明所属类别的标签，还有一个基因集列表。

每次检验都返回一个P值，也可以输出观察数据的检验统计量，用于所有置换的检验统计量和其它一些输出内容。

还有一些用来绘图的函数。

【注】WWtest统计量用于离散的类似正态的分布；KStest统计量用于离散的类似Smirnov的分布（斯米尔诺夫分布）；MDtest统计量用于离散的类似正态的分布；WGSNCA统计量用于连续的类似正态的分布；Fisher-combinedF-统计量用于连续的类似卡方的分布

目前存在许多用于检验差异表达基因集的方法，并且都有共同的名字叫基因集分析（genesetanalysis）。

基因集分析（GSA）可以分为竞争（competitive）分析和自包含（self-contained）分析。

竞争分析方法比较一个基因集与它的补集，而自包含方法在两条件下比较一个基因集是否是差异表达基因集。

一些竞争分析方法受到基因组覆盖度和数据过滤等影响，不相关的数据甚至噪声[4]可以都可能增加检验效能。

两个条件下表达上调和表达下调的基因在基因集中的比例也能影响其它的竞争分析方法[5]。

鉴于以上问题，GSAR包只使用自包含的基因集方法。

通过自包含方法使用不同的检验统计量构造出不同的统计假设的可能性使得可以制定和探索不同的生物学假说[2]。

对于基因集分析（GSA）来说，除了均值表达向量是否相等以外，假设检验仍然被较少研究。

GSAR包提供了一系列检验原假设和特定备择假设的方法，例如差异分布（WWtest函数），位移或者均值的差异（KStest和MDtest函数），尺度或者方差的差异（RKStest，RMDtest和ArrgFtest函数），还有网络相关性结构的差异（GSNCAtest函数）。

可以在GSAR包中使用的大多数统计方法都是基于网络或者图的（除了GSNCAtest和AggrFtest）。

GSAR使用来自igraph包中的具有丰富功能的igraph数据类型来处理图[6]。

GSAR调用一些来自igraph包中的函数来实现自身的一些方法并且使用igraph类型的绘图方法对生成的图进行可视化。

在本手册的举例和案例分析中使用的数据集来自以下数据包：

ALL，GSVAdata和tweeDEseqCountData。

GSAR本身包含了一个经过预处理的数据集来展示分析流程，在基因集网络相关性分析（GSNCA）方法的论文[1]也使用了这个数据集。

运行例子和案例分析的其它必要软件包还有MASS，GSEABase，annotate，org.Hs.eg.db，genefilter，hgu95av2.db，edgeR。

分析从加载GSAR包开始：

>library（GSAR）

接下来我们介绍公式中使用的数学符号的意义。

考虑到两种不同的生物学表型，表示第一种表型中的样本量，表示第二种表型中的样本量。

每个样本是一个-维的向量，是基因的数量（构成了一个单独的基因集）。

因此，第一个表型下的数据是一个的矩阵，第二个表型下的数据是一个的矩阵，行表示基因列表示样本。

两个表型的样本各自由两个随机向量和来表示。

两个向量是独立同分布（iid）并且有分布函数为和，-维向量的均值是和，维的协方差矩阵是和。

2.最小生成树（Minimumspanningtrees）

2.1第一最小生成树（FirstMST）

可以用一个边加权图来表示合并两个条件和下的多元（-维）观测值，是图中结点的集合。

样本网络中的每一个结点对应一个观测（样本）并且是节点对之间的边的集合。

和下的完全图拥有个结点和个边。

边的权重用观测（样本）对的欧氏距离来估计。

在这里，最小生成树（MST）的定义是所有结点之间的边的子集，并且边长度的和最短。

图中的每个节点对应一个来自和下的-维观测。

MST通过对多元观测相对应的结点位置的秩次提供了对多元观测排序的方法。

排序的目的是为了得到观测（样本）的排序差异和它们的距离（）之间的强相关性。

排序算法可以通过专门的设计以确保特定的备择假设的检验效能[2]。

GSAR包中有五种基于最小生成树的MST假设检验方法：

WWtest，KStest，MDtest，RKStest和RMDtest。

下面的例子使用来自MASS包的随机多元正态数据生成器生成了一个具有20个特征和40个观测的特征集，之后用该数据构建了一个图对象并使用igraph包中的函数得到了它的最小生成树。

>library（MASS）

>set.seed（123）

>nf<-20

>nobs<-60

>zero_vetor<-array（0,c（1,nf））

>cov_mtrx<-diag（nf）

>dataset<-mvrnorm（nobs,zero_vetor,cov_mtrx）

>Wmat<-as.matrix（dist（dataset,method=“euclidean”,

+diag=TRUE,upper=TREU,p=2））

>gr<-graph_from_adjacency_matrix（Wmat,weight=TRUE,

+mode=“undirected”）

>MST<-mst（gr）

MST对象的最小生成树图

2.2第二最小生成树（SMT2）用于相关性和PPI网络（SMT2forcorrelationandPPInetworks）

第二最小生成树的定义是图的最小生成树，也就是完全图的边减去第一最小生成树的边的图再次进行最小生成树分析的结果。

在这里我们用MSTs和MST2来分别表示第一最小生成树和第二最小生成树。

如果结点之间所有边都考虑进去的话（完整网络）MST2中的每个结点的度最小为2。

相关性（共表达）网络定义为一个边加权图，图中的结点V是特征的集合每个结点对应一个特征（基因），图中的边E是结点之间的连接集合，边的权值是通过一些相关性距离度量方法得到的估计值。

这里相关性距离的定义是，和分别表示相关性距离和基因和的相关系数[1]。

相关性网络的第二最小生成树（MST2）给出了基因间的必要连接（互作）的最小集合，我们把它解释成一个功能互作的网络。

在基因集中如果一个基因与大多数的其它基因高度相关，那么这个基因趋于占据中心的地位并且在第二最小生成树中有相对较高的度；因为第一和第二最小生成树中连接结点的最短路径都倾向于经过这个基因。

相比之下，一个有较低的基因间相关性的基因更可能在MST2中占据非中心的位置并且度只有2。

MST2的这个性质使得它成为一个有价值的图形可视化工具，通过突出最有影响力的基因来检验由基因表达数据得到的完整的相关性网络。

例如，上个例子中数据集的MST2可以用以下命令来获得：

##TheinputoffindMST2mustbeamatrixwithrowsand

##columnsrespectivelycorrespondingtogenesand

##columns.

##Therefore,datasetmustbetransposedfirst.

>dataset<-aperm（dataset,c（2,1））

>MST2<-findMST2（dataset）

使用findMST2函数获得MST2

蛋白质相互作用网络（PPInetwork）的定义是一个二元图（binarygraph），是图中结点的集合，每个结点对应一个蛋白质；是结点之间的连接的集合，如果和两个蛋白质之间存在互作关系那么，否则。

一个蛋白质相互作用网络的MST2给出了蛋白质之间的必要互作的最小集合。

文献[7]表明函数findMST2.PPI揭示了具有高度连通的蛋白质在簇中占据中心位置的精细网络结构。

3.统计方法（Statisticalmethods）

GSAR包提供了七种统计方法用来检验与原假设相对的五种特定的备择假设（参见图1），两个函数用来寻找相关性网络和蛋白互作网络的MST2，一个封装的函数用来对MST2画图，一个封装函数用来轻松地对一个基因集列表应用特定的统计方法。

3.1瓦尔德-沃尔夫维茨检验（Wald-Wolfowitztest）

Wald-Wolfowitz（WW）检验的原假设是，备择假设是。

当基因数量时，一元WW检验从对来自两个表型的观测（样本）的升序排序开始并且标记每个样本的类标签（表型）。

然后，计算运行次数（），其中是相同标签的联系序列。

检验统计量是运行次数的函数，它是一个渐近正态分布的。

文献[3]中的多元增广（）是基于MST的。

同单变量情况类似，WW检验的多元增广情况下MST图中的不同类标签的两个结点的所有边都被剔除并且计算余下的互斥树（）的数量[3]。

检验统计量是子树数量的标准化值：

通过大规模置换（）观测（样本）类标签得到的零分布（nulldistribution）并且计算出每一次置换的值。

这个零分布是个近似的正态分布。

P-值的计算如下：

这里是次置换的检验统计量，是观测的检验统计量，是指示器函数。

WWtest函数实现这个检验。

3.2柯尔莫哥洛夫-斯米尔诺夫检验（Kolmogorov-Smirnovtests）

当基因数量时，一元KS检验对来自两类表型的观测进行升序排序。

观测值被排序后其定量值的计算为：

这里（）是对（）排序后小于的观测数量，。

检验统计量是来自X和Y的样本秩的累积分布函数（CumulativeDistributionFunction，CDFs）的最大绝对差值，。

通过大规模置换（）观测（样本）类标签得到的零分布（nulldistribution）并且计算出每一次置换的值。

KS统计量近似服从Smirnov分布并且检验单尾假设。

P-值的计算如下：

其中是次置换的检验统计量，是观测的统计量，是一个指示方程。

排序方案可以被设计成服从一个特定备择假设一边获得更高的检验效能。

两个概率情况：

第一，如果零假设是而备择假设是，那么MST根植于最大距离的结点并且其它结点根据预先定义的高秩次来横向遍历整个树[3]，可以使用HDP.ranking函数来进行计算。

函数KStest实现这个特殊的检验。

第二，如果零假设是而备择假设是，那么MST根植于最小距离的结点（中心）并且从根算起具有最大深度的结点具有更高的秩次序。

因此，MST中其它结点的秩次从根结点开始放射状地增加。

一个树中放射状的结点秩可以用函数radial.ranking来计算。

函数RKStest实现这个特殊的检验。

之前例子中的MST网络展示在下图中第一组中的结点标记为绿色第二组中的结点标记为黄色。

图中的结点的秩次使用HDP.ranking和radial.ranking分别对MST的结点进行HDP和放射横向遍历来获得秩次。

>HDP.ranking（MST）

[1]291916544428342060154217723272548404411283222

[25]43102464731185637353955575093335149511591453

[49]1465845381326523021236

>radial.ranking（MST）

[1]65524474235731183756173433494050935232541452

[25]58759602028551392784548141541230211511225344

[49]38324619243131026162936

随机数据的最小生成树

3.3平均偏差检验（Meandeviationtests）

平均偏差（meandeviation，MD）统计量计算样本和秩次的累积分布函数（CDFs）之间的偏差。

一个大小为的基因集的MD统计量定义如下：

其中：

;

是样本在MST中的秩并且指数设置为0.25以赋给秩次一个适当的权重。

通过大规模置换（）观测（样本）类标签得到统计量的零分布并且计算每一次置换的值。

MD统计量是一个渐近正态分布并且检验双尾假设。

P-值的计算如下：

其中是次置换的检验统计量，是观测的检验统计量，并且是一个指示方程。

类似于KS统计量，适当的排序方案与这个统计量组合能够允许进行特定的备择假设的检验，也就是差异均值（均值偏差检验函数MDtest）和方差差异（放射均值偏差方程RMDtest）

3.4累积F-检验（AggregatedF-test）

一元F-检验用来检测来自一个大小为的基因集中的基因个体的方差差异，并且这些单独的P-值，用将费舍尔概率组合（FisherProbabilityCombining）方法累积起来得到这个基因集的打分统计量（）：

显著性的通过置换样本标签来进行估计并且在大规模置换（）次数下计算值。

P-值的计算如下：

其中是次置换的检验统计量，是观测检验统计量，并且是一个指示方程。

3.5基因集网络相关性分析（Genesetnetcorrelationsanalysis）

3.5.1方法

基因集网络相关性分析（GSNCA）是一个两样本非参数多元差异共表达检验（two-samplenonparametricmultivariatedifferentialcoexpressiontest），它考虑的是特征（基因）间的相关性结构。

给定一个大小为的基因集，检验分配权重系数，给一个条件（表型）下的基因并同时调整这些权重以至于每个基因间的权重达到相等并对这些加权的绝对值相关系数（）求和：

，

这个问题被当作一个特征向量问题来解决并且有唯一解基因间的相关系数矩阵的特征向量相应的最大特征值（细节参见[1]）

检验统计量是给定的两个条件下尺度化的权重向量和的第一范式（每个向量都乘以它的范式）：

这个检验统计量检验零假设其备择假设则是。

这个方法的表现在文献[1]中得到了检验。

P-值的计算使用的方法与之前的WW和KS检验的P-值计算方法一样。

尺度变换的额和大概落在区间[0.5,1.5]，在这个基因集中具有较大值的基因是与其它基因高度相关的。

3.5.2零标准偏差的问题（Theproblemofzerostandarddeviation）

在特殊的案例中，在一个或者两个条件（表型）下集合中的一些特征可能是常数或者在各个样本中的值近似常数水平。

这样的情况几乎不会在芯片数据中出现，但是RNA-seq的计数数据中可能出现一个基因在许多样本中都是0个计数，这表示了改基因没有表达。

这就导致了该基因是0标准偏差或者只有很小的标准偏差。

这样的情况下调用cor命令来计算特征间相关系数的时候程序会报错。

。

为了避免这样的情况，提前核查标准偏差（默认行为）并且如果发现任何的标准偏差小于提前设定的阈值（默认是1e-3）的话，运行程序停止、报错并返回出现问题的特征数量（只有一个特征出现问题也会报错）。

为了把GSNCA应用在计数数据上，我们必须从集合中剔除这些出现问题的特征。

在两个条件下如果一个特征在一些样本中有近似常数的情况出现，置换样本标签可以依概率地把这些样本集中在同一个条件下，可以让阈值筛选出这个很小的标准偏差。

为了允许检验跳过这些置换结果而不耽误程序的运行，可以设置一个允许跳过的特征数量的上限（默认是10）。

一旦数量超过上限，程序就会返回错误信息。

如果使用者能够确定被检验的特征集不包含这样的近似0标准偏差的特征（比如经过过滤的芯片数据），则可以跳过对微小标准偏差的核查步骤来减少程序的执行时间。

4.处理RNA-seq计数数据的注释（NotesonhandlingRNA-seqcountdata）

RNA-seq数据由计数的整数构成，这些整数通常用离散的泊松分布（Poisson）或者负二项分布（NegativeBinomial）来表示。

因此，为芯片数据设计的检验方法（芯片数据服从连续的正态分布）不能直接应用到RNA-seq的数据上。

GSAR包中的非参数检验不需要预先的分布假设并且可以应用在经过适当标准化后的RNA-seq的计数数据上。

标准化应该考虑到样本间的差异（测序库的大小或者测序深度）样本内的擦怀疑（组要是基因长度）。

RPKM值（readsperkilobasepermillion）就是这样的标准化。

然而，因为缺少全面的性能的研究，以下两点必须声明：

用来对RNA-seq数据建模的泊松分布和负二项分布的方差是它们均值的函数。

因此，当对RNA-seq数据使用方差检验的时候就需要高度依赖于标准化的方法。

这个特别的问题到现在也没有被较好的研究和描述。

RNA-seq数据集通常都有许多0，因此，在基因集中至少有一个基因的标准偏差是0的问题很常见并且阻止计算GSNCA必须的相关系数。

一个解决方案是丢弃任何可能产生0或者微小的标准偏差的基因然后在剩余的基因上应用GSNCA方法。

5.个案研究（Casestudies）

这一部分将展示GSAR包对基因集的典型应用流程。

案例分析中我们使用了2个芯片数据集和一个RNA-seq数据集。

5.1p53数据集（Thep53dataset）

5.1.1简介（Introduction）

p53是一个主要的肿瘤抑制蛋白。

p53数据集中有50个NCI-60细胞系分化的样本，17个细胞系是TP53野生型（WT），33个细胞系携带了TP53的突变（MUT）[8,9]。

使用hgu95av2芯片平台获得转录谱，数据可以在BroadInstitute网站上获取（http:

//www.broadinstitute.org/gsea/datasets.jsp）。

5.1.2数据过滤和标准化（Filteringandnormalization）

经过预处理的p53数据集的矩阵包含在GSAR包中。

这个p53数据集是从BroadInstitute网站上下载得到的。

探针水平的信号强度用分位数标准化方法处理并且经过对数尺度的转换。

原始探针的每个Affymetrix探针名称和唯一的基因名称相对应。

那些没有对应到entrez和基因名称的探针都被剔除。

考查具有重复的信号强度的探针，我们选择在WT和MUT两个条件下具有T-统计量最大绝对值的那些探针。

最后，矩阵的行用基因名称标识，列则用样本属于WT