ロトカ=ヴォルテラのような時系列変化をしているものの相関を調べたいんだけど、という相談を受けた。
それで、既にGeneNetというパッケージに目をつけているらしいのでこれをやってみる。
他、ggm, GGMselectというのもあるらしいけど時間があれば…たぶんない。
データセットecoliは、ヒト酸化還元酵素を発現させたあとの102遺伝子の発現をマイクロアレイで調べたデータである。
Description
This data set describes the temporal expression of 102 genes of E. Coli after induction of the expression of SOD (recombinant human superoxide dismutase).
時間は0から数えて9つタイムポイントを設定し、データはにしてある。基準はtime 0なので、ここは0の値になっている。
Format
caulobacter is a longitudinal object containing the data from the Schmidt-Heck et al. (2004) experiment. Essentially, this is a matrix with with 102 columns (=genes) and 9 rows (=time points).
All expression levels are given in log2-ratios with respect to the first time point (i.e. the induction at time 0).
longitudinalというデータ型があるみたいだがようわからん。
is.longitudinal(ecoli) as.longitudinal # longitudinal に変換
library(GeneNet) library(gplots) data(ecoli) cols <- greenred(20) par(mar=c(4.5, 3.5, 3, 1)) image(t(ecoli), axes=FALSE, col=cols) axis(1, at=seq(0, 1, length=ncol(ecoli)), labels=colnames(ecoli), las=2, cex.lab=0.3, tick=FALSE) axis(2, at=seq(0, 1, length=nrow(ecoli)), labels=rownames(ecoli), las=2, tick=FALSE) axis(3, at=seq(0, 1, length=ncol(ecoli)), labels=seq(ncol(ecoli)), las=2, cex.lab=0.3, tick=FALSE)
inferred.pcor <- ggm.estimate.pcor(ecoli) test.results <- ggm.test.edges(inferred.pcor) # 有向の場合は direct=TRUE net <- extract.network(test.results, cutoff.ggm=0.9) res0 <- fdrtool(sm2vec(inferred.pcor), statistic="correlation") res0$param
Estimating optimal shrinkage intensity lambda (correlation matrix): 0.1804
Estimate (local) false discovery rates (partial correlations): Step 1... determine cutoff point Step 2... estimate parameters of null distribution and eta0 Step 3... compute p-values and estimate empirical PDF/CDF Step 4... compute q-values and local fdr Step 5... prepare for plotting
net pcor node1 node2 pval qval prob 1 0.23185664 51 53 2.220446e-16 3.612205e-13 1.0000000 2 0.22405545 52 53 2.220446e-16 3.612205e-13 1.0000000 3 0.21507824 51 52 2.220446e-16 3.612205e-13 1.0000000 4 0.17328863 7 93 3.108624e-15 3.792816e-12 0.9999945 5 -0.13418892 29 86 1.120811e-09 1.093997e-06 0.9999516 6 0.12594697 21 72 1.103836e-08 8.978558e-06 0.9998400 7 0.11956105 28 86 5.890921e-08 3.853588e-05 0.9998400 8 -0.11723897 26 80 1.060525e-07 5.816169e-05 0.9998400 ...
res0$param cutoff N.cens eta0 eta0.SE kappa kappa.SE [1,] 0.03553068 4352 0.9474623 0.005656465 2043.377 94.72265
相関、選んだ2つのノード、p値とq値、probability that edge is nonzero (= 1-local fdr)(?)が返ってくる。
有向グラフにすれば有向グラフでのp値、q値も返ってくる。
確率の高いものを抜き出して元データをプロットすると、それなりに相関が高そうな発現データが得られそう。
par(mfrow=c(2, 5)) plot(ecoli, unlist(net[1:5, paste("node", 1:2, sep="")]), lwd=3, autolayout=FALSE)