psandovalv — Jan 16, 2013, 8:25 PM
######################################################################
# This note can be directly run in R.
# Requires GeneNet 1.2.0 (May 2007)
#######################################################################
# This reproduces the four networks displayed in
# Opgen-Rhein and Strimmer (2006a,b)
# Opgen-Rhein, R., and K. Strimmer. 2006a. Using regularized dynamic
# correlation to infer gene dependency networks from time-series
# microarray data. Proceedings of WCSB 2006 (June 12-13, 2006, Tampere, Finland).
# Opgen-Rhein, R., and K. Strimmer. 2006b. Inferring gene dependency
# networks from genomic longitudinal data: a functional data approach.
# REVSTAT 4:53-65.
#load GeneNet library
require(corpcor)
Loading required package: corpcor
Warning: package 'corpcor' was built under R version 2.15.1
require(fdrtool)
Loading required package: fdrtool
Warning: package 'fdrtool' was built under R version 2.15.1
require(longitudinal)
Loading required package: longitudinal
library("GeneNet")
# get T cell data
data(tcell)
tc44 <- combine.longitudinal(tcell.10, tcell.34)
# estimate partial correlations
pc1 <- ggm.estimate.pcor(tc44, lambda=0) # static, no shrinkage
Specified shrinkage intensity lambda (correlation matrix): 0
pc2 <- ggm.estimate.pcor(tc44, method="dynamic", lambda=0) # dynamic, no shrinkage
Specified shrinkage intensity lambda (correlation matrix): 0
pc3 <- ggm.estimate.pcor(tc44) # static, with shrinkage
Estimating optimal shrinkage intensity lambda (correlation matrix): 0.0209
pc4 <- ggm.estimate.pcor(tc44, method="dynamic") # dynamic, with shrinkage
Estimating optimal shrinkage intensity lambda (correlation matrix): 0.032
# significant edges (local fdr, 0.2 cutoff)
# static, no shrinkage
t1.edges <- ggm.test.edges(pc1)
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
t1.net <- extract.network(t1.edges) # prob > 0.8
Significant edges: 6
Corresponding to 0.36 % of possible edges
t1.net
pcor node1 node2 pval qval prob
1 0.4806 11 28 2.125e-07 0.0002589 0.9997
2 0.4744 40 56 3.191e-07 0.0002589 0.9821
3 -0.4141 31 45 1.134e-05 0.0061309 0.8674
4 0.3635 33 48 1.385e-04 0.0493381 0.8674
5 0.3539 23 28 2.133e-04 0.0632622 0.8674
6 0.3444 46 49 3.217e-04 0.0768025 0.8674
# dynamic, no shrinkage
t2.edges <- ggm.test.edges(pc2)
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
t2.net <- extract.network(t2.edges) # prob > 0.8
Significant edges: 9
Corresponding to 0.54 % of possible edges
t2.net
pcor node1 node2 pval qval prob
1 0.5196 11 28 4.548e-09 7.247e-06 0.9821
2 0.3972 31 33 1.490e-05 1.020e-02 0.9821
3 0.3888 40 56 2.325e-05 1.206e-02 0.9821
4 0.3817 18 44 3.365e-05 1.340e-02 0.9779
5 0.3749 23 28 4.754e-05 1.515e-02 0.9318
6 -0.3544 5 26 1.291e-04 2.979e-02 0.9318
7 0.3503 23 32 1.560e-04 3.300e-02 0.9318
8 0.3477 46 49 1.759e-04 3.504e-02 0.9079
9 0.3415 30 40 2.337e-04 4.138e-02 0.8790
# static, with shrinkage
t3.edges <- ggm.test.edges(pc3)
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
t3.net <- extract.network(t3.edges) # prob > 0.8
Significant edges: 10
Corresponding to 0.6 % of possible edges
t3.net
pcor node1 node2 pval qval prob
1 0.4024 40 56 2.762e-07 0.0003308 0.9997
2 0.3968 11 28 4.158e-07 0.0003308 0.9464
3 0.3164 46 49 7.146e-05 0.0276747 0.9464
4 0.3128 23 28 8.711e-05 0.0303085 0.9464
5 0.3094 33 48 1.053e-04 0.0327667 0.9464
6 -0.2961 31 45 2.119e-04 0.0407350 0.9464
7 0.2959 31 33 2.142e-04 0.0408401 0.9464
8 0.2955 6 45 2.192e-04 0.0410630 0.9464
9 0.2940 29 36 2.363e-04 0.0417716 0.9070
10 0.2892 5 30 3.012e-04 0.0473959 0.9070
# dynamic, with shrinkage
t4.edges <- ggm.test.edges(pc4)
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
t4.net <- extract.network(t4.edges) # prob > 0.8
Significant edges: 21
Corresponding to 1.27 % of possible edges
t4.net
pcor node1 node2 pval qval prob
1 0.4450 11 28 1.152e-10 1.765e-07 0.9981
2 0.3426 18 44 1.249e-06 9.571e-04 0.9951
3 0.3217 23 28 5.750e-06 2.581e-03 0.9951
4 0.3177 40 56 7.618e-06 2.918e-03 0.9811
5 0.2987 23 32 2.736e-05 7.484e-03 0.9811
6 0.2924 46 49 4.115e-05 9.383e-03 0.9811
7 0.2911 31 33 4.461e-05 9.764e-03 0.9231
8 0.2785 37 51 9.701e-05 1.848e-02 0.9231
9 0.2671 5 30 1.899e-04 2.940e-02 0.9231
10 -0.2661 32 41 2.013e-04 3.047e-02 0.9231
11 0.2625 4 55 2.463e-04 3.424e-02 0.9231
12 0.2556 46 57 3.628e-04 4.166e-02 0.9231
13 -0.2554 5 26 3.670e-04 4.188e-02 0.9231
14 0.2528 30 40 4.222e-04 4.453e-02 0.9231
15 0.2518 25 42 4.460e-04 4.556e-02 0.8718
16 0.2463 3 16 5.997e-04 5.457e-02 0.8718
17 0.2458 10 16 6.134e-04 5.528e-02 0.8581
18 0.2387 38 58 8.908e-04 6.825e-02 0.8581
19 0.2383 6 45 9.092e-04 6.898e-02 0.8581
20 0.2381 49 56 9.171e-04 6.929e-02 0.8581
21 -0.2368 40 45 9.838e-04 7.178e-02 0.8250
######## produce plots using graphviz ###########
# produce dot file (without edge labels)
node.labels <- colnames(tc44)
ggm.make.dot(filename="net1.dot", t1.net, node.labels, main="Static")
ggm.make.dot(filename="net2.dot", t2.net, node.labels, main="Dynamic")
ggm.make.dot(filename="net3.dot", t3.net, node.labels, main="Static + Shrink")
ggm.make.dot(filename="net4.dot", t4.net, node.labels, main="Dynamic + Shrink")
# call graphviz to produce a nice graph
system("fdp -T svg -o net1.svg net1.dot") # SVG format
system("fdp -T svg -o net2.svg net2.dot") # SVG format
system("fdp -T svg -o net3.svg net3.dot") # SVG format
system("fdp -T svg -o net4.svg net4.dot") # SVG format
######## produce plots using Rgraphviz ###########
node.labels <- colnames(tc44)
gr1 <- ggm.make.graph( t1.net, node.labels, drop.singles=TRUE)
Loading required package: graph
Warning: package 'graph' was built under R version 2.15.2
gr2 <- ggm.make.graph( t2.net, node.labels, drop.singles=TRUE)
gr3 <- ggm.make.graph( t3.net, node.labels, drop.singles=TRUE)
gr4 <- ggm.make.graph( t4.net, node.labels, drop.singles=TRUE)
gr1
A graphNEL graph with undirected edges
Number of Nodes = 11
Number of Edges = 6
gr2
A graphNEL graph with undirected edges
Number of Nodes = 15
Number of Edges = 9
gr3
A graphNEL graph with undirected edges
Number of Nodes = 16
Number of Edges = 10
gr4
A graphNEL graph with undirected edges
Number of Nodes = 30
Number of Edges = 21
# plot networks
library("Rgraphviz")
Warning: package 'Rgraphviz' was built under R version 2.15.1
Loading required package: grid
par(mfrow=c(2,2))
plot(gr1, main="Static, no shrinkage")
plot(gr2, main="Dynamic, no shrinkage")
plot(gr3, main="Static, with shrinkage")
plot(gr4, main="Dynamic, with shrinkage")
par(mfrow=c(1,1))