Rとオミクス解析：オントロジー、ネットワーク解析が出来る統合解析「dnetパッケージ」

source("http://bioconductor.org/biocLite.R")
biocLite(c("hexbin","ape","supraHex","graph","Rgraphviz","igraph","Biobase","limma","survival","foreach","doMC"))
install.packages("dnet",repos="http://cran.r-project.org",type="source")

#ライブラリの読み込み
library(dnet)

#データの読み込みと設定
data(Fang)
data <- Fang
fdata <- as.data.frame(Fang.geneinfo[, 2:3])
rownames(fdata) <- Fang.geneinfo[, 1]
pdata <- as.data.frame(Fang.sampleinfo[, 2:3])
rownames(pdata) <- Fang.sampleinfo[, 1]

#解析に必要な遺伝子発現データオブジェクトの作成
colmatch <- match(rownames(pdata), colnames(data))
rowmatch <- match(rownames(fdata), rownames(data))
data <- data[rowmatch,c olmatch]
eset <- new("ExpressionSet", exprs = as.matrix(data), phenoData = as(pdata,"AnnotatedDataFrame"), featureData = as(fdata," AnnotatedDataFrame"))

#遺伝子発現値の調整
prob2gene <- function(eset){
  fdat <- fData(eset)
  tmp <- as.matrix(unique(fdat[c("EntrezGene", "Symbol")]))
  forder <- tmp[order(as.numeric(tmp[, 1])),]
  forder <- forder[!is.na(forder[, 1]),]
  rownames(forder) <- forder[, 2]
  system.time({
    dat <- exprs(eset)
    edat <- matrix(data = NA, nrow = nrow(forder), ncol = ncol(dat))
    for (i in 1:nrow(forder)){
      ind <- which(fdat$EntrezGene == as.numeric(forder[i, 1]))
      if (length(ind) == 1){
        edat[i,] <- dat[ind,]
      }else{
        edat[i,] <- apply(dat[ind, ], 2, mean)
      }
    }
  })
  
  rownames(edat) <- rownames(forder)
  colnames(edat) <- rownames(pData(eset))
  esetGene <- new('ExpressionSet', exprs = data.frame(edat), phenoData = as(pData(eset), "AnnotatedDataFrame"), featureData = as(data.frame(forder), "AnnotatedDataFrame"))
  return(esetGene)
}
esetGene <- prob2gene(eset)
esetGene

#igraph用のデータ作成
org.Hs.string <- dRDataLoader(RData = 'org.Hs.string')
ind <- match(V(org.Hs.string)$symbol, rownames(esetGene))

#遺伝子発現値データの整形
esetGeneSub <- esetGene[ind[!is.na(ind)],]
esetGeneSub

#ネットワークデータの作成
nodes_mapped <- V(org.Hs.string)$name[!is.na(ind)]
network <- dNetInduce(g = org.Hs.string, nodes_query = nodes_mapped, 
                      knn = 0, remove.loops = T, largest.comp = T)
V(network)$name <- V(network)$symbol

#ネットワーク解析データの準備,処理に時間がかかります。
D <- as.matrix(exprs(esetGene))
D <- D - as.matrix(apply(D,1,mean),ncol=1)[,rep(1,ncol(D))]
fdr <- dSVDsignif(data=D, num.eigen=NULL, pval.eigen=1e-2, signif="fdr", orient.permutation="row", num.permutation=200, fdr.procedure="stepup", verbose=T)
g <- dNetPipeline(g=network, pval=fdr, method="customised", nsize=30)

#ネットワーク図の作成準備
glayout <- layout.fruchterman.reingold(g)

#図書式の設定
com <- spinglass.community(g, spins=25)
com$csize <- sapply(1:length(com),function(x) sum(com$membership==x))
vgroups <- com$membership
colormap <- "yellow-darkorange"
palette.name <- visColormap(colormap=colormap)
mcolors <- palette.name(length(com))
vcolors <- mcolors[vgroups]
com$significance <- dCommSignif(g, com)

#ノードの設定
vdegrees <- igraph::degree(g)
tmp <- data.frame(ind=1:vcount(g), vgroups, vdegrees)
ordering <- tmp[order(vgroups,vdegrees),]$ind

#図のプロット
visNetArc(g, ordering=ordering, labels=V(g)$geneSymbol, vertex.label.color=vcolors, vertex.color=vcolors, vertex.frame.color=vcolors, vertex.size=log(vdegrees)+0.1, vertex.label.cex=0.4)
ordering <- tmp[order(vgroups,vdegrees),]$ind