Nothing
####### INTERNAL FUNCTIONS
networkmatrix = function(network, nodes, nettype) {
if ((nettype %in% c('directed','undirected')) == FALSE) stop("nettype must be either 'directed' or 'undirected' ")
requireNamespace("igraph", quietly = TRUE)
# adjacency matrix: first converted to igraph
if (inherits(network, 'matrix') & nettype == 'directed') {
colnames(network) = nodes
network = igraph::graph.adjacency(network,"directed", diag=F)
}
if (inherits(network, 'matrix') & nettype == 'undirected') {
colnames(network) = nodes
network = igraph::graph.adjacency(network,"undirected", diag=F)
}
# then, get edgelist from igraph
if (inherits(network, 'igraph')) {
if (is.null(igraph::V(network)$name)) igraph::V(network)$name = nodes
output = igraph::get.edgelist(network, names = T)
}
else if (inherits(network, "dgCMatrix")) {
requireNamespace("Matrix", quietly = TRUE)
indexes = as.matrix(Matrix::summary(network))
output = matrix(nrow = dim(indexes)[1], ncol = 2)
v1 = sort(unique(indexes[,1]))
v2 = sort(unique(indexes[,2]))
output[,1] = as.character(factor(indexes[,1], labels=nodes[v1]))
output[,2] = as.character(factor(indexes[,2], labels=nodes[v2]))
}
return(output)
}
# compute pvalue
pvalue = function(x, ss, K, N) { #2*min(p>,p<)+p=
p1 = phyper(x, m = K, n = N-K, k = ss, lower.tail = FALSE) # p > nab
p2 = phyper(x-1, m = K, n = N-K, k = ss, lower.tail = TRUE) # p < nab
p3 = 1-p1-p2
pval = 2*min(p1,p2)+p3
overenr = (p1<p2)
return(list('p'=pval,'overenr'=overenr))
}
########## neat CLASS
neatc = function(x,...) UseMethod('neat', x)
neat = function(alist, blist = NULL, network, nettype,
nodes, alpha = NULL, mtc.type = 'fdr',
anames = NULL, bnames = NULL) {
if (is.null(alpha) == FALSE) {if ( alpha<=0 | alpha>=1 ) stop('alpha must be in (0,1)')}
if ((nettype %in% c('directed','undirected')) == FALSE) stop("nettype must be either 'directed' or 'undirected' ")
if (nettype == 'directed' & is.null(blist)) stop("blist cannot be null when nettype == 'directed' (see manual)")
# first case: igraph object
if (inherits(network, 'igraph')) {
requireNamespace("igraph", quietly = TRUE)
net = networkmatrix(network, nodes, nettype)
}
# second case A: adjacency matrix
else if (inherits(network, "matrix") & ncol(network)>2) {
if (isSymmetric(network)==TRUE & nettype == 'directed') {
warning('The adjacency matrix is symmetric. Should you set nettype = "undirected"?')
}
if (isSymmetric(network)==FALSE & nettype == 'undirected') {
warning('The adjacency matrix is not symmetric. Should you set nettype = "directed"?')
}
net = networkmatrix(network, nodes, nettype)
}
# second case B: sparse adjacency matrix (class "dgCMatrix")
else if (inherits(network, "dgCMatrix")) {
requireNamespace("Matrix", quietly = TRUE)
net = networkmatrix(network, nodes, nettype)
}
# third case: two-column matrix with labels
else if (inherits(network, "matrix") & ncol(network)==2) net=network
# NB: from this point on, 'net' is the network matrix to be used!!!
if (is.factor(nodes) == T) {nodes = as.character(nodes)}
oa = numeric(); ib = numeric(); nab = numeric()
alogic = vector("list", length(alist))
netred = vector("list", length(alist))
p = numeric(); expect=numeric()
o = vector(); concl = vector(); from = vector(); to = vector()
# UNDIRECTED NETWORKS:
if (nettype == 'undirected') {
eina = vector("list", length(alist))
alogic2 = alogic
if ( is.null( names(alist) ) ) names(alist) = anames
for (i in 1:length(alist)) {
if (is.factor(alist[[i]]) == T) {alist[[i]] = as.character(alist[[i]])}
eina[[i]] = (net[,1] %in% alist[[i]])|(net[,2] %in% alist[[i]])
if (sum(eina[[i]]) == 1) netred[[i]] = as.matrix(t(net[eina[[i]],]))
else if (sum(eina[[i]]) == 0) stop(paste('There are no edges connected to genes in',
names(alist)[i],'. The test cannot be computed.'))
else netred[[i]] = net[eina[[i]],]
alogic[[i]] = (netred[[i]][,1] %in% alist[[i]])
alogic2[[i]] = (netred[[i]][,2] %in% alist[[i]])
oa[i] = sum(alogic[[i]])+sum(alogic2[[i]]) # total degree of A (undirected net!)
}
# first case: each A vs each B (blist provided!)
if (is.null(blist) == FALSE) {
blogic = vector("list", length(blist))
blogic2 = blogic
einb = vector("list", length(blist))
if ( is.null( names(blist) ) ) names(blist) = bnames
for (i in 1:length(blist)) {
if (is.factor(blist[[i]]) == T) {blist[[i]] = as.character(blist[[i]])}
einb[[i]] = (net[,1] %in% blist[[i]])|(net[,2] %in% blist[[i]])
# controls added from version 1.1:
if (sum(einb[[i]]) == 0) stop(paste('There are no edges connected to genes in',
names(blist)[i],'. The test cannot be computed.'))
else if (sum(einb[[i]]) == 1) {
temp = as.matrix(t(net[einb[[i]],]))
blogic[[i]] = (temp[,1] %in% blist[[i]])
blogic2[[i]] = (temp[,2] %in% blist[[i]])
}
else {
blogic[[i]] = (net[einb[[i]],][,1] %in% blist[[i]])
blogic2[[i]] = (net[einb[[i]],][,2] %in% blist[[i]])
}
ib[i] = sum(blogic[[i]]) + sum(blogic2[[i]]) # total degree of B (undirected net!)
}
k=1
for (i in 1:length(alist)) {
for (j in 1:length(blist)) {
from[k] = names(alist)[i]
to[k] = names(blist)[j]
link = (alogic[[i]] & netred[[i]][,2] %in% blist[[j]])|(alogic2[[i]] & netred[[i]][,1] %in% blist[[j]])
nab[k] = sum(link)
temp = pvalue(nab[k], ss = oa[i], K = ib[j], N = 2*dim(net)[1])
p[k] = temp$p
expect[k] = round(ib[j] * oa[i] / (2*dim(net)[1]), 4)
o[k] = temp$overenr
if (is.null(alpha) == FALSE) {
if (p[k] > alpha) concl[k] = 'No enrichment'
else if (p[k] <= alpha & o[k]==T) concl[k] = 'Overenrichment'
else if (p[k] <= alpha & o[k]==F) concl[k] = 'Underenrichment'
}
k = k+1
}
}
}
# second case: blist NOT provided
if (is.null(blist) == TRUE) {
ib = oa # undirected: A = B (and no distinction between indegree and outdegree)
k=1
for (i in 1:(length(alist)-1)) {
for (j in (i+1):length(alist)) {
from[k] = names(alist)[i]
to[k] = names(alist)[j]
link = (alogic[[i]] & netred[[i]][,2] %in% alist[[j]])|(alogic2[[i]] & netred[[i]][,1] %in% alist[[j]])
nab[k] = sum(link)
temp = pvalue(nab[k], ss = oa[i], K = ib[j], N = 2*dim(net)[1])
p[k] = temp$p
expect[k] = round(ib[j] * oa[i] / (2*dim(net)[1]), 4)
o[k] = temp$overenr
if (is.null(alpha) == FALSE) {
if (p[k] > alpha) concl[k] = 'No enrichment'
else if (p[k] <= alpha & o[k]==T) concl[k] = 'Overenrichment'
else if (p[k] <= alpha & o[k]==F) concl[k] = 'Underenrichment'
}
k = k+1
}
}
}
} # end for part for undirected networks
# DIRECTED NETWORKS:
if (nettype == 'directed') {
if ( is.null( names(alist) ) ) names(alist) = anames
for (i in 1:length(alist)) {
if (is.factor(alist[[i]]) == T) {alist[[i]] = as.character(alist[[i]])}
alogic[[i]] = (net[,1] %in% alist[[i]])
oa[i] = sum(alogic[[i]]) # outdegree of A
if ( sum(alogic[[i]]) ==1 ) netred[[i]] = as.matrix(t(net[alogic[[i]],]))
else if (sum(alogic[[i]]) == 0) stop(paste('There are no arrows going out from genes in',
names(alist)[i]),'. The test cannot be computed.')
else netred[[i]] = net[alogic[[i]],]
}
# directed: only relevant case is each A vs B (blist provided!)
if (is.null(blist) == FALSE) {
blogic = vector("list", length(blist))
if ( is.null( names(blist) ) ) names(blist) = bnames
for (i in 1:length(blist)) {
if (is.factor(blist[[i]]) == T) {blist[[i]] = as.character(blist[[i]])}
blogic[[i]] = (net[,2] %in% blist[[i]])
ib[i] = sum(blogic[[i]]) # indegree of B
}
k=1
for (i in 1:length(alist)) {
for (j in 1:length(blist)) {
from[k] = names(alist)[i]
to[k] = names(blist)[j]
nab[k] = sum( netred[[i]][,2] %in% blist[[j]])
temp = pvalue(nab[k], ss = oa[i], K = ib[j], N = dim(net)[1])
p[k] = temp$p
expect[k] = round(ib[j] * oa[i] / dim(net)[1], 4)
o[k] = temp$overenr
if (is.null(alpha) == FALSE) {
if (p[k] > alpha) concl[k] = 'No enrichment'
else if (p[k] <= alpha & o[k]==T) concl[k] = 'Overenrichment'
else if (p[k] <= alpha & o[k]==F) concl[k] = 'Underenrichment'
}
k = k+1
}
}
}
} # end of part for directed networks
# final common code
p.adj = p.adjust(p, method = mtc.type)
if (is.null(alpha) == FALSE) {
# concl was based on p until version 1.1.3, here this is
# changed and based on the adjusted p-value
concl[p.adj > alpha] = 'No enrichment'
results = data.frame(from, to, nab, expect, p, p.adj, concl)
names(results) = c('A', 'B', 'nab', 'expected_nab',
'pvalue', 'adjusted_p', 'conclusion')
}
else {
results = data.frame(from, to, nab, expect, p, p.adj)
names(results) = c('A', 'B', 'nab', 'expected_nab',
'pvalue', 'adjusted_p')
}
class(results)=c('neat','data.frame')
results
}
# summary, plot and print:
summary.neat = function(object, ...) {
cat("Number of comparisons:",dim(object)[1],"\n")
cat("Enrichments at 1% level:",sum(object$adjusted_p<0.01),"\n")
cat("Enrichments at 5% level:",sum(object$adjusted_p<0.05),"\n")
cat("Kolmogorov-Smirnov test for uniformity of p-values:",round(ks.test(x=object$pvalue, y='punif')$p.value,4),"\n")
}
plot.neat = function(x, nbreaks=10, ...) {
par(mfrow=c(1,2))
hist(x$pvalue,breaks=nbreaks, xlim=c(0,1), xlab='p-value', main='Histogram of p-values')
plot(ecdf(x$pvalue),main='P-p plot of p-values',xlab='x',ylab ='F(x)',xlim=c(0.037,0.963),ylim=c(0.037,0.963))
abline(v=1);abline(a=0,b=1,col='red')
legend(0.6,0.2,c('Distribution of','p-values','Uniform','distribution'),col=c('black','white','red','white'),lty=1,cex=0.6)
}
print.neat = function(x,nrows=10,...) {
class(x) = 'data.frame'
if (nrows == 'ALL') {
nrows = dim(x)[1]
}
head(x,nrows)
}
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.