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

ゲノミクスやメタボロミクスなどの生化学分野のみならず、ウェブ解析にもオミクス解析が必要ではと感じています。個人的にはオミクス解析は「１つの事柄（反応）のみではなく、事象として表現（現れる）結果を現時点で使用できる技術を最大限に生かし丸ごと解析しよう！そして新しい解釈を導きだそう」を目的とし、ビックデータは「結果を導きだす要因を探そう！データは多いほど正解に近い要因を発見できるよね？もしかしたら、気づいていない要因があるかも？」と考えています。だいぶ乱暴ですがアプローチの仕方がことなると考えています。そんなオミクス解析が出来るdnetパッケージを紹介します。データの準備が大変ですが、参考URLの例をこなすことで使えるようになると思います。

「dnet」パッケージの導入

下記コードを実行することで導入することができます。
参考URL: http://supfam.org/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")

ネットワーク図の描写

オミクス解析で良く描写される、関係性を示すネットワーク図の描写コードを示します。参照URL先には他にも例があるのでぜひご覧ください。データの準備が大変ですが、オミクス解析は非常に価値がある手法だと考えます。ぜひ、dnetパッケージを使用してみてください。

#&#12521;&#12452;&#12502;&#12521;&#12522;&#12398;&#35501;&#12415;&#36796;&#12415;
library(dnet)

#&#12487;&#12540;&#12479;&#12398;&#35501;&#12415;&#36796;&#12415;&#12392;&#35373;&#23450;
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]

#&#35299;&#26512;&#12395;&#24517;&#35201;&#12394;&#36986;&#20253;&#23376;&#30330;&#29694;&#12487;&#12540;&#12479;&#12458;&#12502;&#12472;&#12455;&#12463;&#12488;&#12398;&#20316;&#25104;
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"))

#&#36986;&#20253;&#23376;&#30330;&#29694;&#20516;&#12398;&#35519;&#25972;
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&#29992;&#12398;&#12487;&#12540;&#12479;&#20316;&#25104;
org.Hs.string <- dRDataLoader(RData = 'org.Hs.string')
ind <- match(V(org.Hs.string)$symbol, rownames(esetGene))

#&#36986;&#20253;&#23376;&#30330;&#29694;&#20516;&#12487;&#12540;&#12479;&#12398;&#25972;&#24418;
esetGeneSub <- esetGene[ind[!is.na(ind)],]
esetGeneSub

#&#12493;&#12483;&#12488;&#12527;&#12540;&#12463;&#12487;&#12540;&#12479;&#12398;&#20316;&#25104;
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

#&#12493;&#12483;&#12488;&#12527;&#12540;&#12463;&#35299;&#26512;&#12487;&#12540;&#12479;&#12398;&#28310;&#20633;,&#20966;&#29702;&#12395;&#26178;&#38291;&#12364;&#12363;&#12363;&#12426;&#12414;&#12377;&#12290;
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)

#&#12493;&#12483;&#12488;&#12527;&#12540;&#12463;&#22259;&#12398;&#20316;&#25104;&#28310;&#20633;
glayout <- layout.fruchterman.reingold(g)

#&#22259;&#26360;&#24335;&#12398;&#35373;&#23450;
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)

#&#12494;&#12540;&#12489;&#12398;&#35373;&#23450;
vdegrees <- igraph::degree(g)
tmp <- data.frame(ind=1:vcount(g), vgroups, vdegrees)
ordering <- tmp[order(vgroups,vdegrees),]$ind

#&#22259;&#12398;&#12503;&#12525;&#12483;&#12488;
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

出力例

少しでも、ウェブや実験の解析が楽になりますように！！