ICDS實(shí)現(xiàn)了綜合考慮DNA-methylation, copy number variation, and gene-expression 數(shù)據(jù)，來識(shí)別失調(diào)的癌癥子通路，文章提供了相應(yīng)的R包來實(shí)現(xiàn)文章中的子通路搜索流程
咱先把包安上，來看看怎么用：
install.packages("ICDS")
library(ICDS)
#查看長幫助文檔
browseVignettes("ICDS")
該包主要提供了下方的函數(shù)
-
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
獲取示例數(shù)據(jù)可以使用 getExampleData函數(shù)
library(ICDS)
# 獲取表達(dá)數(shù)據(jù)
exp_data<-GetExampleData("exp_data")
#查看表達(dá)數(shù)據(jù)
exp_data[1:6, 1:6]
#獲取表達(dá)數(shù)據(jù)的樣本標(biāo)簽label 0/1s
# 0 represents the case sample and 1 represents the control sample
label1<-GetExampleData("label1")
#Student’s t-test計(jì)算 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)
# 計(jì)算拷貝數(shù)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值
為了節(jié)省時(shí)間，這里直接用R包中事先跑好的3種p值來進(jìn)行合并操作
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就是你三種類型數(shù)據(jù)整合的結(jié)果，它包含合并后的顯著性，和基因的風(fēng)險(xiǎn)打分
step3
在前一步獲取了每個(gè)基因的風(fēng)險(xiǎn)打分后，就開始利用貪婪搜索算法搜索子通路了，函數(shù)還有其他參數(shù)FindSubPath(zz, Pathway = "kegg", delta = 0.05, seed_p = 0.05,
min.size = 5, out.F = FALSE, out.file = "Subpath.txt")，delta=0.05說明增刪節(jié)點(diǎn)新的通路得分要比原通路打分高5%以上，直到不滿足這個(gè)條件后停止搜索（目的是不讓子通路無限搜索下去）；seed_p是卡種子節(jié)點(diǎn)的閾值，滿足合并p值小于0.05的節(jié)點(diǎn)視為種子節(jié)點(diǎn)；min.size是限定每個(gè)子通路的最少節(jié)點(diǎn)數(shù)，其他參數(shù)來確定是否輸出
#obtain z-score of each gene
require(graphite)
zz<-GetExampleData("zzz")
subpathdata<-FindSubPath(zz)
step4
優(yōu)化子通路：把重疊率過高的子通路合并在在一起，重疊率我們用雅卡爾系數(shù)來表示,也就是代碼中的overlap參數(shù)
subpathdata<-GetExampleData("subpathdata")
keysubpathways<-opt_subpath(subpathdata,zz,overlap=0.6)
step5
擾動(dòng)分析：獲取子通路顯著性，函數(shù)中提及了兩種擾動(dòng)方法，具體可以參看文章,擾動(dòng)次數(shù)最好1000次以上，示例代碼是為了測(cè)試運(yùn)行更快一點(diǎn)，僅僅使用100次
keysubpathways<Permutation(keysubpathways,zz,nperm1=100,method1=TRUE,nperm2=100,method2=FALSE)
step6
給子通路繪圖:
- subpID是結(jié)果中的子通路基因
- pathway.name 通路名
- zz就是先前步驟獲得的每個(gè)基因的風(fēng)險(xiǎn)打分
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)
掃描下方二維碼關(guān)注生信客部落公眾號(hào)：
