1概述
此软件包实现信号通路影响分析(Signaling Pathway Impact Analysis，SPIA)算法。SPIA使用一组差异表达基因及其倍数变化的信息，以及通路拓扑结构来评估这些通路在研究条件下的重要性。出于说明的目的，当前版本的SPIA算法仅包括用于人类和小鼠物种的过时的KEGG信号通路数据。然而，该软件包的当前版本包括从KEGG xml(KGML)文件生成所需的最新处理通路数据的函数，许可用户可以从KEGG的ftp站点为感兴趣的物种下载这些文件。此外，这些文件也可以使用Download KEGML按钮从每个路径的网页上单独下载。对于给定的生物物种，将被处理和分析的路径是：i)包含SPIA所考虑的基因/蛋白质之间的至少一种关系，以及ii)没有反应。
为SPIA分析准备的过时的KEGG数据于2012年9月7日从KEGG的网站下载的，包括人类和小鼠物种，并且做了预处理。
2使用SPIA软件包进行通路分析
我们使用基于Affymetrix基因芯片并可通过GEO获得的结直肠癌数据集(GSE4107)来演示该软件包的功能。实验包括10个正常样本和12个大肠癌样本。用affy软件包对原始数据进行RMA预处理，用limma软件包进行两组适度化t检验。从limma中的topTable函数获得最终结果的数据框作为SPIA的输入数据。这个名为top的数据框包含在SPIA软件包中的colorectalcancer数据集中。
library(SPIA)
data(colorectalcancer)
options(digits=3)
head(top)为了让SPIA发挥作用，我们需要一个矢量，包含在两组之间具有log2倍数变化并且被认为有差异表达的基因。该矢量的名称必须是Entrez基因ID。以下几行将在top数据框中额外添加一列，将每个Affymetrix探针注释到Entrez ID。由于同一Entrez ID可能有多个探测集，因此有两种简单的方法可以获得每个基因的一个对数倍数变化。第一种选择是使用每个基因最显著的样本间的倍数变化，第二种选择是对同一基因的所有探针的对数倍数变化进行平均。在下面的示例中，我们使用了前一种方法。这个例子中的基因被称为差异表达，只要它们的FDR调整后的p值(q值)小于0.05。下列代码从top数据框开始，生成spia函数输入所需的两个矢量：
library(hgu133plus2.db)
x <- hgu133plus2ENTREZID
top$ENTREZ<-unlist(as.list(x[top$ID]))
top<-top[!is.na(top$ENTREZ),]
top<-top[!duplicated(top$ENTREZ),]
tg1<-top[top$adj.P.Val<0.1,]
DE_Colorectal=tg1$logFC
names(DE_Colorectal)<-as.vector(tg1$ENTREZ)
ALL_Colorectal=top$ENTREZDE_Colorecal是一个包含癌症和正常样本之间差异表达基因log2倍数变化的矢量，而ALL_Colorecal是一个在微阵列上描述所有基因的Entrez ID的矢量。DE_Colorecal的名称是与计算的对数倍数变化相对应的Entrez基因ID。
DE_Colorectal[1:10]
ALL_Colorectal[1:10]SPIA算法将上述两个向量作为输入，并生成从最重要到最不重要的通路的表格。这可以通过调用spia函数来实现，如下所示：
# pathway analysis based on combined evidence;
# use nB=2000 or more for more accurate results
res=spia(de=DE_Colorectal,all=ALL_Colorectal,organism="hsa",
nB=2000,plots=FALSE,beta=NULL,combine="fisher",verbose=FALSE)
#make the output fit this screen
res$Name=substr(res$Name,1,10)
#show first 20 pathways, omit KEGG links
res[1:20,-12]如果在上面的函数调用中将plots参数设置为TRUE，则会为存在差异表达基因的每条通路生成如图1所示的图。这些打印保存在当前目录下的pdf文件中。
根据过度表达证据和基于扰动的证据，可以使用函数plotP获得路径重要性的总体图，如图2所示。colorectal cancer pathway路径以绿色显示。
plotP(res,threshold=0.05)
points(I(-log(pPERT))~I(-log(pNDE)),
data=res[res$ID=="05210",],col="green",pch=19,cex=1.5)图2：colorectalcancer数据集的SPIA证据图。每条路径由一个点表示。用Fisher方法结合pPERT和pNDE得到的全局p值pG经Bonferroni校正后，红色斜线右侧的路径是显著的。在对全局p值pG进行FDR校正之后，蓝色斜线右侧的路径是重要的。
在此图中，水平轴表示p值(负对数)，该p值对应于偶然获得给定路径上至少观察到的基因数量(NDE)的概率。垂直轴表示与恰好在给定路径上获得观测到的总累积量(tA)或更极端的概率相对应的p值(负对数)。Tarca等人描述了pPERT的计算。在图2中，每条路径都显示为一个圆点，而在Bonferroni校正之后5%(由plotP中的阈值参数设置)的有效路径以红色显示。
组合pPERT和pNDE的默认方法是Fisher的乘积法，如Tarca等人(2009)所述。
或者，可以使用正态反演方法结合这两种类型的证据，该方法在pPERT和pNDE同时较低时给出较小的pG值。这与Fisher的方法不同，Fisher的方法在两个p值中只有一个是低的时，就会产生较小的pG值。要使用正态反演方法，可以在调用spia函数时设置参数combine="norminv"，或者通过从spia函数生成的结果数据框开始重新计算pG值。下面介绍了后一种方法，其中调用函数combfunc。
res$pG=combfunc(res$pNDE,res$pPERT,combine="norminv")
res$pGFdr=p.adjust(res$pG,"fdr")
res$pGFWER=p.adjust(res$pG,"bonferroni")
plotP(res,threshold=0.05)
points(I(-log(pPERT))~I(-log(pNDE)),data=res[res$ID=="05210",],
col="green",pch=19,cex=1.5)图3：结直肠癌数据集的SPIA证据图。每条路径由一个点表示。用正态反演方法结合pPERT和pNDE得到的全局p值pG经Bonferroni校正后，红色曲线右侧的路径是显著的。在对全局p值pG进行FDR校正之后，蓝色曲线右侧的路径是重要的。
由具有以红色突出显示的差异表达基因的KEGG提供的路径图像可以通过在网络浏览器中粘贴由函数spia产生的数据框的KEGGLINK列中可用的链接来获得。例如,
res[,"KEGGLINK"][20]
#[1] "http://www.genome.jp/dbget-bin/show_pathway?hsa04662+4792+2932+10892+10000+5579+10451+5291+5293+8503+1380+975+2213+2353+3725+3845+4893+25780+4772+4773+4776+11261+5533+5534+63928+5879+5880+974+973+5594+6654+6655+118788"是上面的res数据框中显示第20条路径的KEGG链接。注意，这些数据集的结果可能与Tarca等人描述的结果不同。因为a)这里使用的路径数据库被更新，并且b)默认Beta值被改变。SPIA使用的描述基因/蛋白质之间不同类型关系(如激活或抑制)的图形的有向邻接矩阵可在用于人类的extdata/hsaSPIA.RData文件中获得。SPIA考虑的关系类型和给予它们的默认权重(β系数)为：
rel<-c("activation","compound","binding/association","expression","inhibition",
"activation_phosphorylation","phosphorylation","inhibition_phosphorylation",
"inhibition_dephosphorylation","dissociation","dephosphorylation",
"activation_dephosphorylation","state change","activation_indirect effect",
"inhibition_ubiquination","ubiquination", "expression_indirect effect",
"inhibition_indirect effect","repression","dissociation_phosphorylation",
"indirect effect_phosphorylation","activation_binding/association",
"indirect effect","activation_compound","activation_ubiquination")
beta=c(1,0,0,1,-1,1,0,-1,-1,0,0,1,0,1,-1,0,1,-1,-1,0,0,1,0,1,1)
names(beta)<-rel
cbind(beta)给定关系类型的值为0会导致从所有通路的分析中丢弃这些类型的关系。用户可以通过设置spia函数调用的beta参数随时更改beta的默认值。
用户可以生成自己的基因/蛋白质关系数据，并将其以hsaSPIA.RData文件中所示的列表格式放入。在此文件中，每个路径数据都包含在一个列表中：
load(file=paste(system.file("extdata/hsaSPIA.RData",package="SPIA")))
names(path.info[["05210"]])path.info[["05210"]][["activation"]][25:35,30:40]在上面的矩阵中，只允许0和1值。1表示列给出的基因/蛋白与行给出的基因/蛋白具有“激活”型关系。使用其他R软件包，如graph和Rgraphviz，可以可视化每条途径中每种类型的基因/蛋白质关系的丰度。首先，我们加载所需的包并创建一个函数，该函数可用于将SPIA使用的任何路径的每种类型的关系绘制为图形。
library(graph)
library(Rgraphviz)
plotG<-function(B){
nnms<-NULL;colls<-NULL
mynodes<-colnames(B)
L<-list();
n<-dim(B)[1]
for (i in 1:n){
L[i]<-list(edges=rownames(B)[abs(B[,i])>0])
if(sum(B[,i]!=0)>0){
nnms<-c(nnms,paste(colnames(B)[i],rownames(B)[B[,i]!=0],sep="~"))
}
}
names(L)<-rownames(B)
g<-new("graphNEL",nodes=mynodes,edgeL=L,edgemode="directed")
plot(g)
}然后，基于hsaSPIA数据，我们绘制了ErbB信号通路中的“激活”关系。
plotG(path.info[["04012"]][["activation"]])3用SPIA解析最新的KEGG XML文件
这里我们假设用户从KEGG ftp站点获得给定物种的所有感兴趣通路的KGML(xml)文件(或从KEGG网站逐个下载它们)。作为示例，我们在SPIA包安装的extdata/keggxml/hsa文件夹中包含了四个这样的文件，以演示如何解析这些文件并在路径的结果集合上运行SPIA。
mydir=system.file("extdata/keggxml/hsa",package="SPIA")
dir(mydir)
makeSPIAdata(kgml.path=mydir,organism="hsa",out.path="./")
res<-spia(de=DE_Colorectal, all=ALL_Colorectal, organism="hsa",data.dir="./")
res[,-12]http://www.advaitabio.com/.提供了名为PathwayGuide的SPIA商业版本，其中包括可视化、速度和用户界面方面的附加功能。
