Nothing
# Calculate the Gini index of <x>
Gini <- function(x)
{
n <- length(x)
x <- sort(x)
1/(n-1)*(n+1-2*(sum((n+1-seq_len(n))*x)/sum(x)))
}
# Compare the Hamming distances of the successor states of states and their perturbed copies.
testTransitionRobustness <- function(network, accumulate=TRUE, params=list())
{
res <- do.call("perturbTrajectories", c(list(network=network, measure="hamming"), params))
if (accumulate)
return(res$value)
else
return(res$stat)
}
# Test attractor robustness by searching the original attractor
# in <params$numSamples> perturbed copies of <network> or by perturbing the state trajectories
# and checking whether the attractors change.
testAttractorRobustness <- function(network, accumulate=TRUE, params=list())
{
origAttrs <- getAttractors(network,canonical=TRUE)
# decode parameters
if (is.null(params$perturb))
perturb <- "functions"
else
perturb <- params$perturb
if (params$perturb == "trajectories")
{
params$perturb <- NULL
if (!is.null(params$copies))
{
params$numSamples <- params$copies
params$copies <- NULL
}
res <- do.call("perturbTrajectories", c(list(network=network, measure="attractor"), params))
if (accumulate)
return(res$value)
else
return(res$stat)
}
else
{
if (is.null(params$method))
method <- "bitflip"
else
method <- params$method
if (is.null(params$maxNumBits))
maxNumBits <- 1
else
maxNumBits <- params$maxNumBits
if (is.null(params$numStates))
numStates <- max(1,2^length(network$genes)/100)
else
numStates <- params$numStates
if (is.null(params$simplify))
simplify <- (perturb[1] == "states")
else
simplify <- params$simplify
if (is.null(params$readableFunctions))
readableFunctions <- FALSE
else
readableFunctions <- params$readableFunctions
if (is.null(params$excludeFixed))
excludeFixed <- TRUE
else
excludeFixed <- params$excludeFixed
if (!is.null(params$numSamples))
copies <- params$numSamples
else
if (is.null(params$copies))
copies <- 100
else
copies <- params$copies
perturbationResults <- unlist(sapply(seq_len(copies),function(copy)
{
# get attractors of perturbed network
perturbedAttrs <- getAttractors(perturbNetwork(network,
perturb=perturb,
method=method,
maxNumBits=maxNumBits,
numStates=numStates,
simplify=simplify,
readableFunctions=readableFunctions,
excludeFixed=excludeFixed))
# try to find original attractors in perturbed network
attractorIndices <- sapply(origAttrs$attractors,function(attractor1)
{
index <- which(sapply(perturbedAttrs$attractors,function(attractor2)
{
identical(attractor1,attractor2)
}))
if (length(index) == 0)
NA
else
index
})
return(sum(!is.na(attractorIndices)))
}))
if (accumulate)
# return overall percentage of found attractors
return(sum(perturbationResults)/(length(origAttrs$attractors) * copies) * 100)
else
# return percentage of found attractors for each run
return(perturbationResults/(length(origAttrs$attractors)) * 100)
}
}
# Calculate the in-degrees of states in the network.
# If <accumulate> is true, the in-degrees are accumulated
# using the Gini coefficient
testIndegree <- function(network,accumulate=TRUE,params)
{
attr <- getAttractors(network)
graph <- plotStateGraph(attr,plotIt=FALSE)
# calculate in-degree using igraph
if (accumulate)
# accumulate using Gini index
return(Gini(degree(graph,mode="in",loops=TRUE)))
else
# return the raw degrees
return(degree(graph,mode="in",loops=TRUE))
}
# Calculate the Kullback-Leibler distance of
# two distributions <x> and <y>.
# <bins> is the number of bins used for discretization.
# <minVal> is the minimum value to be used instead of zero
kullbackLeiblerDistance <- function(x,y,bins=list(),minVal=0.00001)
{
x <- cut(x,breaks=bins,include.lowest=T,right=F)
y <- cut(y,breaks=bins,include.lowest=T,right=F)
tx <- table(x)
ty <- table(y)
tx <- tx/length(x)
ty <- ty/length(y)
tx[tx < minVal] <- minVal
ty[ty < minVal] <- minVal
return(sum(tx*(log(tx/ty))))
}
# Generic function to test properties of <network> against random networks.
# <numRandomNets> specifies the number of random networks to generate.
# <testFunction> is the name of a function that returns a distribution of test values or a test statistic for each network,
# which receives <testFunctionParams> as optional parameters.
# If <accumulation> is "characteristic", the test function must return a characteristic/test statistic value for each value, and
# a histogram of these value is plotted with a line for the original network.
# If <accumulation> is "kullback_leibler", a histogram of the Kullback-Leibler distances of the test value distribution for the original network
# and the random networks is plotted.
# <sign.level> is the desired significance level for <network> in comparison to the random networks.
# <functionGeneration>,<simplify>,<noIrrelevantGenes>,<validationFunction>,<failureInterations,
# <d_lattice>,<zeroBias> are the corresponding parameters of generateRandomNKNetwork().
# <title> is the title of the plot, <xlab> is its x axis caption, <breaks> is the corresponding histogram parameter, and ... supplies further
# graphical parameters
testNetworkProperties <- function(network, numRandomNets=100, testFunction="testIndegree",
testFunctionParams=list(),
accumulation=c("characteristic","kullback_leibler"),
alternative=c("greater","less"),
sign.level=0.05, drawSignificanceLevel=TRUE,
klBins,klMinVal=0.00001,
linkage=c("uniform","lattice"),
functionGeneration=c("uniform","biased"),
validationFunction, failureIterations=10000,
simplify=FALSE, noIrrelevantGenes=TRUE,
d_lattice=1, zeroBias=0.5,
title="", xlab, xlim, breaks=30,
...)
{
stopifnot(inherits(network,"BooleanNetwork") || inherits(network,"SymbolicBooleanNetwork"))
if (is.character(testFunction))
testFunctionName <- testFunction
else
testFunctionName <- ""
testFunction <- match.fun(testFunction)
accumulate <- (match.arg(accumulation) == "characteristic")
origResult <- testFunction(network,accumulate,testFunctionParams)
numGenes <- length(network$interactions)
if (inherits(network,"SymbolicBooleanNetwork"))
inputGenes <- sapply(network$interactions,function(interaction)length(getInputs(interaction)))
else
inputGenes <- sapply(network$interactions,function(interaction)length(interaction$input))
if (missing(validationFunction))
validationFunction <- NULL
randomResults <- lapply(seq_len(numRandomNets),function(i)
{
randomNet <- generateRandomNKNetwork(n=numGenes,
k=inputGenes,
topology="fixed",
linkage=linkage,
functionGeneration=functionGeneration,
validationFunction=validationFunction,
failureIterations=failureIterations,
simplify=simplify,
noIrrelevantGenes=noIrrelevantGenes,
d_lattice=d_lattice,
zeroBias=zeroBias)
randomRes <- testFunction(randomNet,accumulate,testFunctionParams)
return(randomRes)
})
if (accumulate)
randomResults <- unlist(randomResults)
args <- list(...)
res <- switch(match.arg(accumulation,c("characteristic","kullback_leibler")),
characteristic = {
# get one value for each random network, and plot a histogram of these values
if (missing(xlab))
{
xlab <- switch(testFunctionName,
testIndegree = "Gini index of state in-degrees",
testAttractorRobustness = "% of identical attractors",
testTransitionRobustness = "Normalized Hamming distance",
"accumulated results"
)
}
if (missing(xlim))
{
#xlim <- switch(testFunctionName,
# testIndegree = c(0,1),
# testAttractorRobustness = c(0,100),
# c(min(c(origResult,randomResults)),
# max(c(origResult,randomResults)))
#)
xlim <- range(c(origResult, randomResults))
}
alternative <- match.arg(alternative, c("greater","less"))
# calculate p-value
if (alternative == "greater")
pval <- sum(randomResults < origResult) / length(randomResults)
else
pval <- sum(randomResults > origResult) / length(randomResults)
# plot histogram
if (testFunctionName == "testIndegree" | testFunctionName == "testAttractorRobustness")
{
r <- hist(randomResults,xlim=xlim,xlab=xlab,main=title,xaxt="n",...)
axis(side=1,at=seq(xlim[1],xlim[2],length.out=11))
}
else
{
# plot with default axis
r <- hist(randomResults,xlim=xlim,xlab=xlab,main=title,...)
}
# plot result for original network
abline(v=origResult,col="red")
if (alternative == "greater")
text(x=origResult,pos=2,y=max(r$counts)*0.75,
labels=paste("> ",round(pval * 100),"%\nof random results",sep=""),
col="red",cex=0.75)
else
text(x=origResult,pos=4,y=max(r$counts)*0.75,
labels=paste("< ",round(pval * 100),"%\nof random results",sep=""),
col="red",cex=0.75)
if (drawSignificanceLevel)
# plot line for significance level
{
if (alternative == "greater")
{
quant <- quantile(randomResults,1.0-sign.level)
abline(v=quant,col="blue")
text(x=quant,pos=2,y=max(r$counts)*0.85,
labels=paste((1.0-sign.level) * 100,"% quantile",sep=""),
col="blue",cex=0.75)
}
else
{
quant <- quantile(randomResults,sign.level)
abline(v=quant,col="blue")
text(x=quant,pos=4,y=max(r$counts)*0.85,
labels=paste(sign.level * 100,"% quantile",sep=""),
col="blue",cex=0.75)
}
}
list(hist=r,pval=1.0-pval,significant=(1.0-pval<=sign.level))
},
kullback_leibler =
{
# a distribution of values is returned for each network,
# plot the Kullback-Leibler distances to the original network
if (missing(xlab))
xlab <- "Kullback-Leibler distance"
if (missing(klBins))
{
bins <- unique(c(origResult,unlist(randomResults)))
bins <- c(bins,max(bins) + 1)
}
else
{
bins <- unique(c(origResult,unlist(randomResults)))
if (klBins < length(bins))
bins <- seq(min(bins),max(bins),length.out=klBins+1)
else
bins <- c(bins,max(bins) + 1)
}
vals <- sapply(randomResults,function(results)
kullbackLeiblerDistance(origResult,results,bins=bins,minVal=klMinVal))
r <- hist(vals,xlab=xlab,main=title,breaks=breaks,...)
list(hist=r)
},
stop("'accumulation' must be one of \"characteristic\",\"kullback_leibler\""))
return(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.