[R]R包】R包ICDS使用简介：Identification

前阵子在简书上写了【ICDS:】Identification of cancer dysfunctional subpathway by integrating DNA-methylation, copy number variation, and gene-expression data

ICDS实现了综合考虑DNA-methylation, copy number variation, and gene-expression 数据，来识别失调的癌症子通路，文章提供了相应的R包来实现文章中的子通路搜索流程

咱先把包安上，来看看怎么用：

install.packages("ICDS")
library(ICDS)
#查看长帮助文档
browseVignettes("ICDS")

该包主要提供了下方的函数

• getExpp,getMethp,getCnvp ：calculate p-values or corrected p-values for each gene.
• coverp2zscore,combinep_two,combinep_three :Convert p-values or corrected p-values to z-scores.
• FindSubPath :search for interested subpathways in each entire pathway.
• opt_subpathOptimize interested subpathways.
• Permutation :calculate statistical significance for these interested subpathways .
• PlotSubpathway绘制子通路图

step1

获取示例数据可以使用 getExampleData函数

library(ICDS)
# 获取表达数据
exp_data<-GetExampleData("exp_data")
#查看表达数据
exp_data[1:6, 1:6]
#获取表达数据的样本标签label  0/1s
# 0 represents the case sample and 1 represents the control sample
label1<-GetExampleData("label1")

#Student’s t-test计算 p-value label的0 代表 normal samples;1 代表 cancer samples.
exp.p<-getExpp(exp_data,label = label1,p.adjust = FALSE)
label2<-GetExampleData("label2")
meth_data<-GetExampleData("meth_data")
meth.p<-getMethp(meth_data,label = label2,p.adjust = FALSE)

# 计算拷贝数p-value
#obtain Copy number variation data
cnv_data<-GetExampleData("cnv_data")
#obtion amplified genes
amp_gene<-GetExampleData("amp_gene")
#obtion deleted genes
del_gene<-GetExampleData(("del_gene"))
#calculate p-values or corrected p-values for each gene
cnv.p<-getCnvp(exp_data,cnv_data,amp_gene,del_gene,p.adjust=FALSE,method="fdr")
#p.adjust=TURE代表矫正
#p.adjust=FALSE代表不矫正

---

step2

合并三种类型的p值

为了节省时间，这里直接用R包中事先跑好的3种p值来进行合并操作

exp.p<-GetExampleData("exp.p")
meth.p<-GetExampleData("meth.p")
cnv.p<-GetExampleData("cnv.p")
#calculate z-scores for p-values of each kind of data
zexp<-coverp2zscore(exp.p)
zmeth<-coverp2zscore(meth.p)
zcnv<-coverp2zscore(cnv.p)
#combine two kinds of p-values,then,calculate z-score for them
zz<-combinep_two(exp.p,meth.p)
#combine three kinds of p-values,then,calculate z-score for them
zzz<-combinep_three(exp.p,meth.p,cnv.p)
#zzz就是你三种类型数据整合的结果，它包含合并后的显著性，和基因的风险打分

---

step3

在前一步获取了每个基因的风险打分后，就开始利用贪婪搜索算法搜索子通路了，函数还有其他参数FindSubPath(zz, Pathway = "kegg", delta = 0.05, seed_p = 0.05,
min.size = 5, out.F = FALSE, out.file = "Subpath.txt")，delta=0.05说明增删节点新的通路得分要比原通路打分高5%以上，直到不满足这个条件后停止搜索（目的是不让子通路无限搜索下去）；seed_p是卡种子节点的阈值，满足合并p值小于0.05的节点视为种子节点；min.size是限定每个子通路的最少节点数，其他参数来确定是否输出

#obtain z-score of each gene
require(graphite)
zz<-GetExampleData("zzz")
subpathdata<-FindSubPath(zz) 

---

step4

优化子通路：把重叠率过高的子通路合并在在一起，重叠率我们用雅卡尔系数来表示,也就是代码中的overlap参数

subpathdata<-GetExampleData("subpathdata")
keysubpathways<-opt_subpath(subpathdata,zz,overlap=0.6)

---

step5

扰动分析：获取子通路显著性，函数中提及了两种扰动方法，具体可以参看文章,扰动次数最好1000次以上，示例代码是为了测试运行更快一点，仅仅使用100次

keysubpathways<Permutation(keysubpathways,zz,nperm1=100,method1=TRUE,nperm2=100,method2=FALSE)

---

step6

给子通路绘图:

• subpID是结果中的子通路基因
• pathway.name 通路名
• zz就是先前步骤获得的每个基因的风险打分

require(graphite)
require(org.Hs.eg.db)
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)

扫描下方二维码关注生信客部落公众号：

生信客部落