時系列データの相関

MikuHatsune2013-04-09

ロトカ=ヴォルテラのような時系列変化をしているものの相関を調べたいんだけど、という相談を受けた。
それで、既に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つタイムポイントを設定し、データはlog_2にしてある。基準は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)