inst/doc/arabidopsis-net.R

In GeneNet: Modeling and Inferring Gene Networks

# /*
# This is an R script containing R markdown comments.  It can be run as is in R.
# To generate a document containing the formatted R code, R output and markdown 
# click the "Compile Notebook" button in R Studio, or run the command
# rmarkdown::render() - see http://rmarkdown.rstudio.com/r_notebook_format.html
# */

#' ---
#' title: "Arabidopsis Thaliana Network"
#' output: pdf_document
#' author: ""
#' date: Example for GeneNet 1.2.13 (August 2015) or later
#' ---

#' This note reproduces the “Arabidopsis thaliana” network example from
#' R. Opgen-Rhein and K. Strimmer. 2007. *From correlation to causation 
#' networks: a simple approximate learning algorithm and its application 
#' to high-dimensional plant gene expression data.*
#' BMC Syst. Biol. **1**: 37.
#' (http://dx.doi.org/10.1186/1752-0509-1-37)

#'  
#'The original source of the data is 
#' Smith et al. 2004. *Diurnal changes in the transcriptom encoding 
#' enzymes of starch metabolism provide evidence for both transcriptional
#' and posttranscriptional regulation of starch metabolism in Arabidopsis 
#' leaves.*  Plant Physiol. **136**: 2687-2699.
#'  

#' This example was suggested by
#' Papapit Ingkasuwan, Division of Biotechnology,
#' School of Bioresources and Technology,
#' King Mongkut's University of Technology Thonburi, Bangkok, Thailand.


#'
#' # Inspect Data

library("GeneNet")
data("arth800")
summary(arth800.expr)

#' Plot time series:
plot(arth800.expr, 1:9)

#' Inspect pairwise scatter plots:
panel.cor = function(x, y, digits=2, prefix="", cex.cor)
{
    usr = par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r = abs(cor(x, y))
    txt = format(c(r, 0.123456789), digits=digits)[1]
    txt = paste(prefix, txt, sep="")
    if(missing(cex.cor)) cex = 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex * r)
}
#+ fig.width=8, fig.height=10
pairs(arth800.expr[,1:9], lower.panel=panel.smooth, upper.panel=panel.cor)


#'
#' # Compute Partial Correlations and Select Relevant Edges

pcor.dyn = ggm.estimate.pcor(arth800.expr, method = "dynamic")
arth.edges = network.test.edges(pcor.dyn,direct=TRUE)
dim(arth.edges)

#' We use the strongest 150 edges:
arth.net = extract.network(arth.edges, method.ggm="number", cutoff.ggm=150)


#'
#' # Construct Graph

library("graph")

node.labels = as.character(1:ncol(arth800.expr))
gr = network.make.graph( arth.net, node.labels, drop.singles=TRUE) 

#' Some information about the graph

#' Number of nodes:
num.nodes(gr)

#' Correlations:
edge.info(gr)$weight

#' Number of directed ("forward") and undirected ("none") edges:
table(  edge.info(gr)$dir )
# /*
# forward    none 
#     55      95
# */

#'
#' # Well-Connected Nodes

#' Nodes connected with many edges:
sort(node.degree(gr), decreasing=TRUE)[1:10]
# /*
# 570  81 783  47 422 558 452 539 738 272
# 20  17  10   9   9   9   8   8   8   7
# */

arth800.descr[570]
# /* [1] "AP2 transcription factor - like protein" */

arth800.descr[81]
# /* [1] "ATRPAC43; DNA binding / DNA-directed RNA polymerase;  */

arth800.descr[558]
# /* [1]"structural constituent of ribosome; */

arth800.descr[539]
# /* [1] "DNA binding / transcription factor;  */

arth800.descr[783]
# /* [1] "RNA binding / RNA methyltransferase; */



#'
#' # Plot Network

library("Rgraphviz")

#' For a more beautiful plot of the network set node and edge parameters:

#' Set global node and edge attributes:
globalAttrs = list()
globalAttrs$edge = list(color = "black", lty = "solid", lwd = 1, arrowsize=1)
globalAttrs$node = list(fillcolor = gray(.95), shape = "ellipse", fixedsize = FALSE)

#' Set attributes of some particular nodes:
nodeAttrs = list()
nodeAttrs$fillcolor = c('570' = "red", "81" = "red") # highlight hub nodes

#' Set edge attributes:
edi = edge.info(gr) # edge directions and correlations
edgeAttrs = list()
edgeAttrs$dir =  edi$dir # set edge directions 
cutoff = quantile(abs(edi$weight), c(0.2, 0.8)) # thresholds for line width / coloring
edgeAttrs$lty = ifelse(edi$weight < 0, "dotted", "solid") # negative correlation
edgeAttrs$color = ifelse( abs(edi$weight <= cutoff[1]), "grey", "black") # lower 20% quantile
edgeAttrs$lwd = ifelse(abs(edi$weight >= cutoff[2]), 2, 1) # upper 20% quantile

#+ fig.width=8, fig.height=7
plot(gr, attrs = globalAttrs, nodeAttrs = nodeAttrs, edgeAttrs = edgeAttrs, "fdp")