识别肿瘤功能失调子通路的方法ICDS

导语

GUIDE ╲

子通路是指具有特定生物学功能的生物通路的局部区域。随着大规模测序数据的产生使我们有更多的机会来研究癌症发生的分子机制。研究DNA甲基化、拷贝数变异(CNV)和基因表达改变对致瘤的失调子通路分子状态的潜在影响是很必要的。本工作提出一个通过整合多组学数据和通路拓扑信息来识别癌症功能失调子通路(ICDS)的方法。利用肝癌(LIHC)、头颈部鳞状细胞癌(HNSC)、宫颈鳞状细胞癌和宫颈腺癌的数据集，验证了ICDS在识别异常子通路方面的有效性。进一步将ICDS和其他识别子通路的方法）（只考虑DNA甲基化、CNV或基因表达）进行比较，通过这些分析，证实ICDS比其他三种只考虑一种数据类型的方法更能识别癌症相关的子通路。

方法介绍

1.计算Gene Risk Score

首先通过整合以下三种类型的数据来计算gene risk评分：DNA甲基化、CNV和基因表达。

（1）用Student’s t-test 计算差异基因表达、差异甲基化，用GISTIC2方法识别样本的拷贝数扩增和缺失，将样本根据基因的拷贝数状态分组，然后用Student’s t-test 计算差异表达。然后校正P值。

（2）通过使用Fisher’s combined probability test计算上述3种P值的联合P值，来计算gene risk score (RS)

（3）根据反正态累积分布函数(CDF)将P值转化为z-score 。

2.计算Subpathway-Activity Score

对于一个KEGG的通路，使用贪婪搜索算法在判别得分（discriminative scores）局部最大的通路中识别关键失调子通路。

（1）搜索算法从一个具有显著高风险得分(p < 0.001)的种子基因i开始，迭代扩展，然后选择该种子基因的一个邻居来形成当前的子通路。对于一个子通路k, 活性评分(activity score ，ASk)为子通路中成员基因RS的平均值，由下式计算

在每次迭代中，算法都从当前子通路的相邻基因中选取一个基因，使得ASk+1和ASk之间的增长达到最大。当ASk+1评分超过（1+r）ASk时没有其他基因增加，或者在当前子通路中任何两个节点之间的距离大于3，以保持局部搜索，搜索算法将停止。选择的改进率（improvement rate ）r避免了太大的subpathway区域（导致增加了冗余弱信息）。在应用于生物网络的贪心启发式算法中，有证据表明参数r = 0.05是合适的 (Chuang et al., 2007)。当同一通路中每对子通路之间的Jaccard index大于0.6，将这两个子通路结合，这样确保了识别的子通路包含了更多的信息，减少了冗余。此外，只考虑了5个以上基因和不到100个基因的子通路，以避免过窄或过宽的功能子通路。

3. Subpathway的显著性检验

使用置换检验来计算这些关键的失调子通路的统计显著性水平。对每个候选子通路提供了两种统计检验方法：whole gene-based perturbation和local-gene perturbation。用户可以选择测试方法。

（1）第一个测试方法扰动了通路网络中整个基因列表上的基因标签，并重新计算子通路的活性评分，记为ASk_perm1，用于检验真实子通路与疾病表型的相关性。

（2）第二个测试方法扰乱了该子通路所属通路中的基因名称，并重新计算了该子通路的活性评分，记为ASk_perm2，用于检验真实子通路与通路结构的相关性。

使用Benjamini和Hochberg方法校正P值。FDR<0.001作为显著subpathway的阈值。

R包介绍

1.数据进行T检验

所用示例数据：

exp_data：TCGA的样本的表达数据

meth_data：TCGA的样本的甲基化数据

cnv_data：TCGA的样本的拷贝数数据

amp_gene：所识别的扩增基因列表

del_gene：所识别的缺失基因列表

（1）getCnvp

对CNV数据进行t-test

exp_data<-GetExampleData("exp_data")
meth_data<-GetExampleData("meth_data")
cnv_data<-GetExampleData("cnv_data")
amp_gene<-GetExampleData("amp_gene")
del_gene<-GetExampleData(("del_gene"))

getCnvp(exp_data,cnv_data,amp_gene,del_gene,
p.adjust=FALSE,method="fdr")
#若p.adjust=TRUE，则返回校正后的P值

（2）getExpp

对基因表达数据进行t-test

profile<-GetExampleData("exp_data")
label<-GetExampleData("label1")

getExpp(profile,label,p.adjust=FALSE)
#若p.adjust=TRUE，则返回校正后的P值

（3）getMethp

对甲基化数据进行t-test

profile<-GetExampleData("meth_data")
label<-GetExampleData("label2")

getMethp(profile,label,p.adjust=FALSE)
#若p.adjust=TRUE，则返回校正后的P值

2.计算联合P值

使用Fisher’s combined probability test方法计算联合P值

（1）combinep_three

计算联合三个维度的P值，用以计算z-score

所用示例文件（其实就是上述1.（1）（2）（3）中得到的结果文件）：

exp.p：Student’s t-test 计算差异基因表达，然后校正P值

meth.p：Student’s t-test 差异甲基化，然后校正P值

cnv.p：GISTIC2方法识别样本的拷贝数扩增和缺失，将样本根据基因的拷贝数状态分组，然后用Student’s t-test 计算差异表达，然后校正P值

exp.p<-GetExampleData("exp.p")
meth.p<-GetExampleData("meth.p")
cnv.p<-GetExampleData("cnv.p")

combinep_three(exp.p,meth.p,cnv.p)

（2）combinep_two

计算联合两个维度的P值，用以计算z-score

exp.p<-GetExampleData("exp.p")
meth.p<-GetExampleData("meth.p")

combinep_two(exp.p,meth.p)

3.coverp2zscore

对P值计算z-score

exp.p<-GetExampleData("exp.p")
meth.p<-GetExampleData("meth.p")
cnv.p<-GetExampleData("cnv.p")

coverp2zscore(exp.p)
coverp2zscore(meth.p)
coverp2zscore(cnv.p)

4.FindSubPath

使用贪婪搜索算法在每个完整的通路中来检索key subpathways

所用示例数据：zzz，是基因表达P值、校正后P值、z-score数据

require(graphite) #载入graphite包
#graphite包，将通路拓扑转换为基因网络
zz<-GetExampleData("zzz")
k<-FindSubPath(zz,Pathway = "kegg",delta = 0.05, seed_p = 0.05,
               min.size = 5, out.F = TRUE)
#delta搜索子路径每一步的扩散系数
#seed_p将p值小于seed_p的基因定义为种子基因
#min.size子通路的最小尺寸
#out.F是否输出子通路

好像是网络的问题，小编下载不了数据，这里就不展示结果了

5.opt_subpath

opt_subpath对感兴趣的子路径进行优化。如果两个通路共享的基因数量大于每个通路基因的重叠比例，则合并两个通路。

示例文件：subpathdata

zz<-GetExampleData("zzz")
subpathdata<-GetExampleData("subpathdata")
keysubpathways<-opt_subpath(subpathdata,zz,overlap=0.6)
#overlap，每个通路基因的重叠率

输出:

6.Permutation

采用置换检验 的method 1和method 2计算这些最优子通路的统计显著性水平。

require(graphite)
keysubpathways<-GetExampleData("keysubpathways")  
zzz<-GetExampleData("zzz")
keysubpathways<-Permutation(keysubpathways,zzz[,3],nperm1=10,
method1=TRUE,nperm2=10,method2=FALSE)
#nperm1,执行method1的置换次数
#method1,是否使用method1
#nperm2，执行method2的置换次数
#method2，是否使用method2
head(keysubpathways)

7. PlotSubpathway

当用户输入一个基因列表时，绘制一个网络图

require(graphite)
subpID<-unlist(strsplit("G6PC/HK3/GPI/FBP1/ALDOA/G6PC2","/"))
pathway.name="Glycolysis / Gluconeogenesis"
zzz<- GetExampleData("zzz")
PlotSubpathway(subpID=subpID,pathway.name=pathway.name,
               zz=zzz,Pathway="kegg",
               layout = layout.fruchterman.reingold)
#subpID，一个感兴趣的子通路的基因列表
#pathway.name，感兴趣的子通路的名称
#zz，每个基因的z-score
#Pathway，通路数据库的名称
#layout,对布局规范函数的调用

小编总结

在我们的研究中，可能往往无法识别一个与研究工作相关的通路，那么也可能是通路中的某个子通路在发挥作用，而一个通路中的子通路往往更能解释疾病。今天向大家介绍了一个挖掘子通路的工具包，并且同时考虑了基因组、转录组和表观组的影响。所以结合到你的相关工作实践看看吧！

引用：

Liu S, Zheng B, Sheng Y, et al. Identification of Cancer Dysfunctional Subpathways by Integrating DNA Methylation, Copy Number Variation, and Gene-Expression Data. Front Genet. 2019;10:441. Published 2019 May 15. doi:10.3389/fgene.2019.00441