BACKGROUND
nea employs network topology to evaluate functional impact of experimentally determined genes and gene sets by detecting enrichment of previously characterized gene sets, such as pathways.
nea became a natural extension of the gene set enrichment analysis, gsea  <cit>  into the network domain  <cit> .
performance and applicability of gsea have been limited by the following: 1) only a minority of genes possesses specific pathway annotations, 2) it can be applied only to genes altered in a specific way, detectable by the given platform  and not to those regulated via other mechanisms, and 3) statistical power of the analysis is limited by gene set size.
nea largely overcomes these limitations due to one key difference: while gsea counts the number of genes shared between an experimental list and a pathway, nea considers network edges between any genes of the two sets in the global network.
on the other hand, the earlier nea versions had an own drawback: the error rate was estimated in multiple, time-consuming instances of randomized networks.
being a non-parametric approach, the randomization allows solving a wide range of higher-order topological problems.
however, for many applications this would be prohibitively slow.
in the present r package, we implement a much faster parametric estimation of connectivity rate expected by chance and show that this also eliminates bias on small gene sets.
finally, we provide an illustration of how nea increases robustness of non-replicated analyses of gene differential expression.
RESULTS
software functionality
nearender possesses both core and ancillary functions for network enrichment analysis.
it is implemented as an r package and described in details, beyond the official manual, in an r vignette https://cran.r-project.org/web/packages/nearender/vignettes/nearender_vignette.pdf
the input shall contain three components: 1) one or multiple user-defined  gene sets which have to be functionally characterized, called altered gene sets ; 2) a collection of functional gene sets  which would enable functional characterization through their known functions, and 3) a global network of functional coupling .
input can be provided as either text files or pre-processed r lists and matrices.
figure  <dig> presents the relationship between the major components and steps of nearender.fig.
<dig> analysis flow in nearender.
the original matrix of 'omics'  data described a limited number of samples  with a much larger number of gene feature rows.
at the first, preparatory step each sample was described via a characteristic sample-specific altered gene set .
in parallel, a collection of functional gene sets  that share certain functionally annotations  was downloaded or prepared otherwise.
a global gene/protein network  was also provided .
in the course of network enrichment analysis  each ags received as many nea scores as there were fgss, i.e.
obtained coordinates in the multidimensional fgs space.
this created an output matrix of the same number of sample columns but many fewer rows
in fact, although we denote input in terms of “genes”, a range of functional nodes of a biological network can be analyzed with this algorithm, such as protein molecules, genomic regions that encode proteins, micrornas, promoters, and enhancers etc.
nodes listed in ags, fgs, and net should employ the same id format.
the package performs network enrichment analysis through the fast binomial test as described in methods and quantifies enrichment in each ags-fgs pair with a number of statistics: chi-squared score, z-score, p-value, q-value , and respective number of ags-fgs edges expected by chance, calculated with the binomial formula.
in order to run enrichment analysis, the user should prepare the input components.
since the ags is the most dynamic and user-specific part of the input, functionality for ags compilation and processing is most developed.
first, agss can be prepared in advance externally, e.g.
as a list of genomic variations reported in a given genome.
alternatively, the package functions can create agss from sample columns of an r matrix.
this can be done with a number of different algorithms  or pre-defined by user .
in addition, the algorithm ‘toprandom’ generates random agss of a user-set size.
finally, a special function allows direct creation of agss as full sample-specific sets of mutated genes.
fgss are usually imported from a file by listing all members of each functional set.
due to the high network density and, hence, statistical power of nea, there is a special option, which is not available in gsea: single genes can be treated as fgs as well.
a full list of such single-gene fgss can be created from all network nodes of net with parameter as_genes_fgs, so that each fgs item in the output list contains just one gene.
it is practical to use a large pre-compiled fgs collection, such as all ontology terms, pathways, or a union of resources, e.g.
msigdb database  <cit> .
alternatively, users can create custom single- or multi-gene fgs collections of their own.
we note however that using relatively few fgss and/or agss in one analysis  would not allow estimating the q-values properly.
while higher order topological biases are, as discussed above, of arguable importance for nea, another network feature is vital for the enrichment evaluation used in this package.
namely, the parametric algorithm would produce unbiased estimates only in scale-free networks, i.e.
where node connectivity values follow the power law distribution  <cit> .
being an almost ubiquitous feature in the full scale biological networks, this is still not the case in networks that are artificially constructed from e.g.
chip-seq based collections of transcription factor binding events  <cit>  or from computationally predicted microrna-transcript targeting data  <cit> .
if such network components are desirable, we recommend employing software that involves network permutation tests, i.e.
the randomization, such as in  <cit> .
therefore, the package enables evaluating network topology for scale-freeness and second-order dependencies  as well as benchmarking alternative nets using either standard or custom fgss, as described in  <cit> .
briefly, the benchmark consists of as many test cases as there are fgs members in total .
for each such gene, the procedure tests the null hypothesis of the gene not being an fgs member.
the true positive or false negative result is assigned if the gene receives an nea score above or below a certain threshold, respectively.
in parallel, randomly picked genes with close node degree values are tested against the same fgs to estimate specificity via the false positive versus true negative ratio.
the counts of alternative test outcomes tp, tn, fp, and fn at variable nea thresholds are used for plotting roc curves.
comparison to randomization based algorithm
we compared the p-value distributions between the two methods .
the both were capable of adjusting p-values for multiple testing, but we omit the adjusted value analysis here since our gs test set of  <dig> fgs  and chi-squared binomial formula  .
since csb employed a deterministic procedure, we assume that all dispersion of p-value points in the nrz-csb space should be attributed to sampling errors by the stochastic nrz algorithm.
csb and nrz p-values converged sufficiently well only at n >  <dig>  according to the quantile analysis with q-q plots , csb p-values were more conservative than those of nrz.
next, the latter was somewhat less sensitive than csb in regard of small gs, especially of those where both gss were small, below  <dig> or  <dig> nodes altogether .
again, nrz p-values from the computationally feasible but insufficient n = {3;10; 30} exhibited more disagreement with respective csb values.
in addition to the sums n
ags + n
fgs, we stratified the q-q plots by the sum of node degrees c
ags + c
fgs and by the actual number of edges between two gss, n
edges .
the sum c
ags + c
fgs appeared a strongly biasing factor .
the influence of n
edges was almost non-existent.fig.
<dig> comparative sensitivity and sources of bias in randomization-based versus binomial calculation of network enrichment.
network enrichment between all vs. all  <dig> gene sets was analyzed with both nrz and csb methods.
p-value distributions were compared using q-q plots  and scatter plots of log  values .
q-q plots in column  <dig> display both the total distributions , i.e.
regardless of gs size, and distribution fractions that correspond to smaller gs .
the qq-plots in column  <dig> are insets of the black, un-stratified q-q plots of column  <dig>  identity lines  are plotted in dotted grey and dotted red in the qq-plots and scatter plots, respectively.
analogous plots in regard of other factors biasing nrz p-values  are provided in additional file 1: figure s <dig>
since there was no gold standard for network enrichment, we could not unambiguously conclude that csb was more precise than nrz.
however it seems to be preferable because of its, in general, more conservative estimates and, in particular, higher sensitivity on small gs as well as the convergence of nrz at higher values of n. the nrz-specific bias on small gss likely occurred because between-gs connectivity values could not take on negative values, hence the left tail of the distribution of n
edges was compressed and standard deviation must be underestimated.
thus while using nrz, the deviation from csb values and bias due to small gs size were considerably large at n <  <dig>  however, performing many randomizations might make the procedure prohibitively slow.
while using nrz, the computation time grows linearly with the number of randomizations.
for example on a server running  <dig> .32- <dig> el <dig> x86_ <dig> ,one randomization of the given network took around 3 min plus 20 min for counting connectivity in the given fgs collection, so that the total time for  <dig> randomizations was 100* = 2300 min.
for the most direct comparison, the same perl software  <cit>  in the csb mode, i.e.
without network randomization, run 15 min.
in the both cases, the perl process occupied between  <dig> and 400 mb ram.
running r package nearender required less than 5 min and below 200 mb ram.
we note that in many practical tasks the number of either ags or fgs could be significantly reduced.
however as stated above, the current package is mostly meant for generating nea scores that could be used as predictive features in machine learning - and those agsxfgs matrices are supposed to be large.
robustness of results obtained from different replicates
statistically underpowered experimental designs are not welcomed in the scientific community but, due to the lack of resources, are still practiced.
in particular, it is common to represent patients in large cohort screens with single samples  <cit> .
both the gsea and nea analyses measure enrichment in order to summarize signals from individual genes' to the level of pathways and biological processes.
this feature suggests a potential increase in robustness of conclusions in experiments that lack replicates.
we decided to investigate this robustness of using single-gene expression values versus gsea and nea under different scenarios using cell transcriptome data from the fantom <dig> cage rna sequencing set  <cit> .
in order to make results comparable between nea and gsea, we use the binomial version of the latter  <cit> , which allowed applying the analyses to the same input gene sets.
commonly, statistical power is analyzed via extrapolating variance estimates  <cit> .
since gsea and nea do not provide such estimates, we instead watched consistency of conclusions drawn from individual replicates.
the samples, in  <dig> to  <dig> biological  replicates, originated from different fibroblast, epithelial, and smooth muscle cells.
in total,  <dig> replicates from  <dig> cell types quantified rna transcript tags mapped to  <dig> genes.
this provided a broad range of degrees of dissimilarity between type-specific transcriptomes.
the dissimilarity was evaluated as a fraction of significantly de genes in fully replicated designs and ranged from  <dig> %  to  <dig> %  de genes.
this set-up allowed us to model situations of analyzing de of each gene g using single samples, e.g.
g
ai vs. g
bj on samples i and j from cell types a and b, and then using all available replicates for a and b : ga1…na vs. gb1…nb.
note that calculating de p-values was thus possible only in the latter case.
otherwise, de could only be ranked by values of fold change between a
i and b
j. we quantified de in all the  <dig> possible a
i - b
j pairs between the  <dig> cell types.
as an example, if there were n
a =  <dig> and n
b =  <dig> replicates, then we could model de measurements on 5*3 =  <dig> pairs of single samples plus one analysis ga1…na vs. gb1…nb.
in the following, we review preservation of de estimates across different cell samples contrasts analyzed with three different methods: the default de analysis on individual genes, gsea, and nea.
the preservation was evaluated in the form of spearman rank correlation coefficients between: 1) fold change values of individual genes for gene-wise analysis“as is”,2) p-values of fgs enrichment scores obtained with binomial gsea, and 3) p-values of fgs enrichment scores obtained with csb nea.
the example scatterplots at fig.
<dig>  display agreement between the same replicate pairs of the same contrast between fibroblast gingivial vs. tenocyte cells.
the fold change values and p-values, although being different in their nature and scale, allowed us to rank the results in the same way as a biological researcher would prioritize them.
we note here that in gsea and nea neither p-values nor enrichment scores  explicitly provide sampling errors similar to between-sample variance in the de analysis.
all such comparisons between sample pairs were quantified by spearman rank correlations and used as individual points in the plots on the right .
one observation from the plots c, and e was that regardless of the correlation strength, nea possessed a much higher statistical power to detect enrichment – which has been discussed in details by alexeyenko et al.
in  <dig>  <cit> .
indeed, no gsea p-values appeared significant after the bonferroni adjustment for multiple testing in either replicate pair .
for comparison, almost 1/3rd of the nea scores were significant in the both replicate pairs.fig.
<dig> agreement between biological replicates in alternative approaches to differential expression analysis.
a, c, e: examples of spearman rank correlations between de analyses without replicates on gingivial epithelial versus tenocyte cells using samples from two different donor pairs .
since p-values in the non-replicated de analyses were not available, the plot a  represents raw fold change values for  <dig> genes, while plots c and e represent p-values for  <dig> fgss.
in c and e, the vertical and horizontal brown dotted lines delineate p-values significant after the bonferroni correction.
dashed grey line: the linear regression fit of y on x .
nea and gsea values were obtained for agss representing de genes in each analysis that satisfied the criterion abs )> <dig>  i.e.
the 4-fold change in either direction.
these example rank r values from a, c, and e are plotted at respectively b, d, and f as big orange dots.
thus in b, d, and f, the spearman coefficients from a, c, and e  as well as from all other pairwise comparisons are plotted as a function of relative difference strength between respective transcriptomes .
as an example, according to the results of fully powered de analysis for plots a, c, and e,  <dig> % genes were de with adjusted p-value <  <dig> .
hence the value  <dig>  is used as x-coordinate for the orange dots.
the grey dotted linear regression line and the spearman rank r value quantify the relations between x and y. black and colored points correspond to 1) correlation values of non-replicated de analyses with each other and 2) correlations where one of the two analyses was replicated, respectively.
the boxplots g and h summarize results across all cell types and ags versions for non-replicated vs. half-replicated options  and .
all pairwise contrasts of gsea and raw against nea were significant with p <  <dig> .
in additional file 1: figure s <dig>  six plots for nea and gsea represent results of using all the six alternative ags versions
gsea and nea were each tested using six alternative variants of ags, generated for each contrast: a) lists of top  <dig>   <dig>  and  <dig> de genes ranked by fold change and b) lists of genes with absolute log
2
values exceeding  <dig>   <dig>  and  <dig>  ags sizes of the form  were not fixed and depended on the magnitude of the difference between the two transcriptomes.
when using all available replicates in the analyses , the genes were ranked by p-value with functions lmfit and ebayes in r package limma  <cit> .
in  we required that, in addition to the fold change cutoff, the adjusted p-values did not exceed  <dig> .
overall, the scores from nea agreed with each other much better than those from gsea .
by investigating the ags- and method-specific scatterplots for agss of type”abs ) > 2” , we observed that the level of correlation between the analyses  depended on strength of the pairwise difference between cell type transcriptomes, expressed as the fraction of de genes .
ideally, we would observe independence of the difference strength  and perfect correlation between de results from both individual sample pairs and replicated analyses.
however, it was unrealistic to expect large and efficient agss from few  significantly de genes.
still, we can conclude that nea output was considerably closer to the desirable pattern than those of the gene-wise and gsea analyses.
indeed, in addition to the stronger overall correlation, nea performed particularly well on weakly different cell types pairs .
the only exception could be found for agss top <dig> , where gsea appeared superior over nea.
however, this option would be practically unusable in gsea because of its low statistical power on small gene sets and, hence, few significant enrichment values.
we also could compare consistency of results obtained via i) individual sample pairs against each other  and ii) results from individual sample pairs against the replicated analysis .
one can also see that the option , i.e.
the comparison of non-replicated to replicated analyses , often exhibited poorer correlation than comparisons  , especially for gsea .
discussion
we demonstrated that estimates of enrichment from the network randomization procedure were biased.
they could sufficiently converge to respective csb values only at very large, impractical numbers of randomizations runs.
while considering the “lack-of-replicates” scenario, we assumed that the former shall always be preferred.
however, our analysis demonstrated the higher  robustness of nea scores in the full absence or shortage of replicates.
this conclusion about superiority of nea over gene-wise analyses is, of course, only relevant when the research problem can be approached with an exploratory analysis at a pathway level and not when it requires identifying individual consistently deregulated genes.
in addition, this advantage of nea explains the earlier observed higher robustness of its scores, compared to individual gene profiles, as descriptors of experimental models of cancer-inhibitory fibroblasts  <cit> .
importantly, biological replicates of patients’ samples are rarely available in clinical cohorts, so that we anticipate efficient usage of nea in such setups.
nea can also be performed using the online tool available at https://www.evinet.org.
the latter is more interactive, user-friendly, and focused on visualization, while the package nearender enables high performance, possibly with parallel computation.
in combination with the ancillary functions for input preparation and benchmarking nearender is meant to become a practically useful part of bioinformatics and biostatistics r pipelines.
CONCLUSIONS
output of nearender allows using sample- or patient-specific pathway scores as predictive features of phenotype or clinical variables.
compared to gsea, this raises statistical power and robustness to practically acceptable levels.
compared to previous implementations of nea, the analysis runs much faster and removes bias due to smaller gene set size.
package nearender also contains previously not available ancillary functions for preparing and benchmarking input.

 nea employs network topology to evaluate functional impact of experimentally determined genes and gene sets by detecting enrichment of previously characterized gene sets, such as pathways. nea became a natural extension of the gene set enrichment analysis, gsea into the network domain. performance and applicability of gsea have been limited by the following: 1) only a minority of genes possesses specific pathway annotations, 2) it can be applied only to genes altered in a specific way, detectable by the given platform and not to those regulated via other mechanisms, and 3) statistical power of the analysis is limited by gene set size. nea largely overcomes these limitations due to one key difference: while gsea counts the number of genes shared between an experimental list and a pathway, nea considers network edges between any genes of the two sets in the global network. on the other hand, the earlier nea versions had an own drawback: the error rate was estimated in multiple, time-consuming instances of randomized networks. being a non-parametric approach, the randomization allows solving a wide range of higher-order topological problems. however, for many applications this would be prohibitively slow. in the present r package, we implement a much faster parametric estimation of connectivity rate expected by chance and show that this also eliminates bias on small gene sets. finally, we provide an illustration of how nea increases robustness of non-replicated analyses of gene differential expression.  software functionality nearender possesses both core and ancillary functions for network enrichment analysis. it is implemented as an r package and described in details, beyond the official manual, in an r vignette https://cran.r-project.org/web/packages/nearender/vignettes/nearender_vignette.pdf the input shall contain three components: 1) one or multiple user-defined gene sets which have to be functionally characterized, called altered gene sets; 2) a collection of functional gene sets which would enable functional characterization through their known functions, and 3) a global network of functional coupling. input can be provided as either text files or pre-processed r lists and matrices. figure  presents the relationship between the major components and steps of nearender.fig. analysis flow in nearender. the original matrix of'omics' data described a limited number of samples with a much larger number of gene feature rows. at the first, preparatory step each sample was described via a characteristic sample-specific altered gene set. in parallel, a collection of functional gene sets that share certain functionally annotations was downloaded or prepared otherwise. a global gene/protein network was also provided. in the course of network enrichment analysis each ags received as many nea scores as there were fgss, i.e. obtained coordinates in the multidimensional fgs space. this created an output matrix of the same number of sample columns but many fewer rows in fact, although we denote input in terms of “genes”, a range of functional nodes of a biological network can be analyzed with this algorithm, such as protein molecules, genomic regions that encode proteins, micrornas, promoters, and enhancers etc. nodes listed in ags, fgs, and net should employ the same id format. the package performs network enrichment analysis through the fast binomial test as described in methods and quantifies enrichment in each ags-fgs pair with a number of statistics: chi-squared score, z-score, p-value, q-value, and respective number of ags-fgs edges expected by chance, calculated with the binomial formula. in order to run enrichment analysis, the user should prepare the input components. since the ags is the most dynamic and user-specific part of the input, functionality for ags compilation and processing is most developed. first, agss can be prepared in advance externally, e.g. as a list of genomic variations reported in a given genome. alternatively, the package functions can create agss from sample columns of an r matrix. this can be done with a number of different algorithms or pre-defined by user. in addition, the algorithm ‘toprandom’ generates random agss of a user-set size. finally, a special function allows direct creation of agss as full sample-specific sets of mutated genes. fgss are usually imported from a file by listing all members of each functional set. due to the high network density and, hence, statistical power of nea, there is a special option, which is not available in gsea: single genes can be treated as fgs as well. a full list of such single-gene fgss can be created from all network nodes of net with parameter as_genes_fgs, so that each fgs item in the output list contains just one gene. it is practical to use a large pre-compiled fgs collection, such as all ontology terms, pathways, or a union of resources, e.g. msigdb database. alternatively, users can create custom single- or multi-gene fgs collections of their own. we note however that using relatively few fgss and/or agss in one analysis would not allow estimating the q-values properly. while higher order topological biases are, as discussed above, of arguable importance for nea, another network feature is vital for the enrichment evaluation used in this package. namely, the parametric algorithm would produce unbiased estimates only in scale-free networks, i.e. where node connectivity values follow the power law distribution. being an almost ubiquitous feature in the full scale biological networks, this is still not the case in networks that are artificially constructed from e.g. chip-seq based collections of transcription factor binding events or from computationally predicted microrna-transcript targeting data. if such network components are desirable, we recommend employing software that involves network permutation tests, i.e. the randomization, such as in. therefore, the package enables evaluating network topology for scale-freeness and second-order dependencies as well as benchmarking alternative nets using either standard or custom fgss, as described in. briefly, the benchmark consists of as many test cases as there are fgs members in total. for each such gene, the procedure tests the null hypothesis of the gene not being an fgs member. the true positive or false negative result is assigned if the gene receives an nea score above or below a certain threshold, respectively. in parallel, randomly picked genes with close node degree values are tested against the same fgs to estimate specificity via the false positive versus true negative ratio. the counts of alternative test outcomes tp, tn, fp, and fn at variable nea thresholds are used for plotting roc curves. comparison to randomization based algorithm we compared the p-value distributions between the two methods. the both were capable of adjusting p-values for multiple testing, but we omit the adjusted value analysis here since our gs test set of fgs and chi-squared binomial formula. since csb employed a deterministic procedure, we assume that all dispersion of p-value points in the nrz-csb space should be attributed to sampling errors by the stochastic nrz algorithm. csb and nrz p-values converged sufficiently well only at n >  according to the quantile analysis with q-q plots, csb p-values were more conservative than those of nrz. next, the latter was somewhat less sensitive than csb in regard of small gs, especially of those where both gss were small, below or nodes altogether. again, nrz p-values from the computationally feasible but insufficient n = {3;10; 30} exhibited more disagreement with respective csb values. in addition to the sums n ags + n fgs, we stratified the q-q plots by the sum of node degrees c ags + c fgs and by the actual number of edges between two gss, n edges. the sum c ags + c fgs appeared a strongly biasing factor. the influence of n edges was almost non-existent.fig. comparative sensitivity and sources of bias in randomization-based versus binomial calculation of network enrichment. network enrichment between all vs. all gene sets was analyzed with both nrz and csb methods. p-value distributions were compared using q-q plots and scatter plots of log values. q-q plots in column display both the total distributions, i.e. regardless of gs size, and distribution fractions that correspond to smaller gs. the qq-plots in column are insets of the black, un-stratified q-q plots of column identity lines are plotted in dotted grey and dotted red in the qq-plots and scatter plots, respectively. analogous plots in regard of other factors biasing nrz p-values are provided in additional file 1: figure s since there was no gold standard for network enrichment, we could not unambiguously conclude that csb was more precise than nrz. however it seems to be preferable because of its, in general, more conservative estimates and, in particular, higher sensitivity on small gs as well as the convergence of nrz at higher values of n. the nrz-specific bias on small gss likely occurred because between-gs connectivity values could not take on negative values, hence the left tail of the distribution of n edges was compressed and standard deviation must be underestimated. thus while using nrz, the deviation from csb values and bias due to small gs size were considerably large at n <  however, performing many randomizations might make the procedure prohibitively slow. while using nrz, the computation time grows linearly with the number of randomizations. for example on a server running.32- el x86_,one randomization of the given network took around 3 min plus 20 min for counting connectivity in the given fgs collection, so that the total time for randomizations was 100* = 2300 min. for the most direct comparison, the same perl software in the csb mode, i.e. without network randomization, run 15 min. in the both cases, the perl process occupied between and 400 mb ram. running r package nearender required less than 5 min and below 200 mb ram. we note that in many practical tasks the number of either ags or fgs could be significantly reduced. however as stated above, the current package is mostly meant for generating nea scores that could be used as predictive features in machine learning - and those agsxfgs matrices are supposed to be large. robustness of results obtained from different replicates statistically underpowered experimental designs are not welcomed in the scientific community but, due to the lack of resources, are still practiced. in particular, it is common to represent patients in large cohort screens with single samples. both the gsea and nea analyses measure enrichment in order to summarize signals from individual genes' to the level of pathways and biological processes. this feature suggests a potential increase in robustness of conclusions in experiments that lack replicates. we decided to investigate this robustness of using single-gene expression values versus gsea and nea under different scenarios using cell transcriptome data from the fantom cage rna sequencing set. in order to make results comparable between nea and gsea, we use the binomial version of the latter, which allowed applying the analyses to the same input gene sets. commonly, statistical power is analyzed via extrapolating variance estimates. since gsea and nea do not provide such estimates, we instead watched consistency of conclusions drawn from individual replicates. the samples, in to biological replicates, originated from different fibroblast, epithelial, and smooth muscle cells. in total, replicates from cell types quantified rna transcript tags mapped to genes. this provided a broad range of degrees of dissimilarity between type-specific transcriptomes. the dissimilarity was evaluated as a fraction of significantly de genes in fully replicated designs and ranged from % to % de genes. this set-up allowed us to model situations of analyzing de of each gene g using single samples, e.g. g ai vs. g bj on samples i and j from cell types a and b, and then using all available replicates for a and b: ga1…na vs. gb1…nb. note that calculating de p-values was thus possible only in the latter case. otherwise, de could only be ranked by values of fold change between a i and b j. we quantified de in all the possible a i - b j pairs between the cell types. as an example, if there were n a =  and n b =  replicates, then we could model de measurements on 5*3 =  pairs of single samples plus one analysis ga1…na vs. gb1…nb. in the following, we review preservation of de estimates across different cell samples contrasts analyzed with three different methods: the default de analysis on individual genes, gsea, and nea. the preservation was evaluated in the form of spearman rank correlation coefficients between: 1) fold change values of individual genes for gene-wise analysis“as is”,2) p-values of fgs enrichment scores obtained with binomial gsea, and 3) p-values of fgs enrichment scores obtained with csb nea. the example scatterplots at fig.  display agreement between the same replicate pairs of the same contrast between fibroblast gingivial vs. tenocyte cells. the fold change values and p-values, although being different in their nature and scale, allowed us to rank the results in the same way as a biological researcher would prioritize them. we note here that in gsea and nea neither p-values nor enrichment scores explicitly provide sampling errors similar to between-sample variance in the de analysis. all such comparisons between sample pairs were quantified by spearman rank correlations and used as individual points in the plots on the right. one observation from the plots c, and e was that regardless of the correlation strength, nea possessed a much higher statistical power to detect enrichment – which has been discussed in details by alexeyenko et al. in. indeed, no gsea p-values appeared significant after the bonferroni adjustment for multiple testing in either replicate pair. for comparison, almost 1/3rd of the nea scores were significant in the both replicate pairs.fig. agreement between biological replicates in alternative approaches to differential expression analysis. a, c, e: examples of spearman rank correlations between de analyses without replicates on gingivial epithelial versus tenocyte cells using samples from two different donor pairs. since p-values in the non-replicated de analyses were not available, the plot a represents raw fold change values for genes, while plots c and e represent p-values for fgss. in c and e, the vertical and horizontal brown dotted lines delineate p-values significant after the bonferroni correction. dashed grey line: the linear regression fit of y on x. nea and gsea values were obtained for agss representing de genes in each analysis that satisfied the criterion abs )> i.e. the 4-fold change in either direction. these example rank r values from a, c, and e are plotted at respectively b, d, and f as big orange dots. thus in b, d, and f, the spearman coefficients from a, c, and e as well as from all other pairwise comparisons are plotted as a function of relative difference strength between respective transcriptomes. as an example, according to the results of fully powered de analysis for plots a, c, and e, % genes were de with adjusted p-value < . hence the value is used as x-coordinate for the orange dots. the grey dotted linear regression line and the spearman rank r value quantify the relations between x and y. black and colored points correspond to 1) correlation values of non-replicated de analyses with each other and 2) correlations where one of the two analyses was replicated, respectively. the boxplots g and h summarize results across all cell types and ags versions for non-replicated vs. half-replicated options and. all pairwise contrasts of gsea and raw against nea were significant with p < . in additional file 1: figure s six plots for nea and gsea represent results of using all the six alternative ags versions gsea and nea were each tested using six alternative variants of ags, generated for each contrast: a) lists of top and de genes ranked by fold change and b) lists of genes with absolute log 2 values exceeding and ags sizes of the form were not fixed and depended on the magnitude of the difference between the two transcriptomes. when using all available replicates in the analyses, the genes were ranked by p-value with functions lmfit and ebayes in r package limma. in we required that, in addition to the fold change cutoff, the adjusted p-values did not exceed. overall, the scores from nea agreed with each other much better than those from gsea. by investigating the ags- and method-specific scatterplots for agss of type”abs ) > 2”, we observed that the level of correlation between the analyses depended on strength of the pairwise difference between cell type transcriptomes, expressed as the fraction of de genes. ideally, we would observe independence of the difference strength and perfect correlation between de results from both individual sample pairs and replicated analyses. however, it was unrealistic to expect large and efficient agss from few significantly de genes. still, we can conclude that nea output was considerably closer to the desirable pattern than those of the gene-wise and gsea analyses. indeed, in addition to the stronger overall correlation, nea performed particularly well on weakly different cell types pairs. the only exception could be found for agss top, where gsea appeared superior over nea. however, this option would be practically unusable in gsea because of its low statistical power on small gene sets and, hence, few significant enrichment values. we also could compare consistency of results obtained via i) individual sample pairs against each other and ii) results from individual sample pairs against the replicated analysis. one can also see that the option, i.e. the comparison of non-replicated to replicated analyses, often exhibited poorer correlation than comparisons, especially for gsea. discussion we demonstrated that estimates of enrichment from the network randomization procedure were biased. they could sufficiently converge to respective csb values only at very large, impractical numbers of randomizations runs. while considering the “lack-of-replicates” scenario, we assumed that the former shall always be preferred. however, our analysis demonstrated the higher robustness of nea scores in the full absence or shortage of replicates. this conclusion about superiority of nea over gene-wise analyses is, of course, only relevant when the research problem can be approached with an exploratory analysis at a pathway level and not when it requires identifying individual consistently deregulated genes. in addition, this advantage of nea explains the earlier observed higher robustness of its scores, compared to individual gene profiles, as descriptors of experimental models of cancer-inhibitory fibroblasts. importantly, biological replicates of patients’ samples are rarely available in clinical cohorts, so that we anticipate efficient usage of nea in such setups. nea can also be performed using the online tool available at https://www.evinet.org. the latter is more interactive, user-friendly, and focused on visualization, while the package nearender enables high performance, possibly with parallel computation. in combination with the ancillary functions for input preparation and benchmarking nearender is meant to become a practically useful part of bioinformatics and biostatistics r pipelines.  output of nearender allows using sample- or patient-specific pathway scores as predictive features of phenotype or clinical variables. compared to gsea, this raises statistical power and robustness to practically acceptable levels. compared to previous implementations of nea, the analysis runs much faster and removes bias due to smaller gene set size. package nearender also contains previously not available ancillary functions for preparing and benchmarking input. 
however, previous implementations of nea involved multiple runs of network randomization and were therefore slow. that requires, firstly, increasing statistical power and biological relevance via a network enrichment analysis and, secondly, a fast and convenient procedure for rendering the original data into a space of pathway scores. further, we compare estimates from these two alternative procedures and demonstrate that the summarization of individual genes to pathways increases the statistical power compared to both the default differential expression analysis on individual genes and the state-of-the-art gene set enrichment analysis. here, we present a new r package nearender which can transform raw'omics' features of experimental or clinical samples into matrices describing the same samples with many fewer nea-based pathway scores. the package nearender is complementary to the online nea tool evinet and, unlike of the latter, enables high performance of computations off-line. the r package nearender version is available at cran repository however, pathway scores can also be used as predictive features in machine learning. this is done via a parametric estimation of the null binomial distribution and is thus much faster and less biased than randomization procedures. of gene profiles' confluence to the pathway level, allows exploring molecular landscapes using functionally annotated gene sets. the package also contains functions for preparing input, modeling null distributions, and evaluating alternative versions of the global network.
 the statistical evaluation of pathway enrichment, i.e. of gene profiles' confluence to the pathway level, allows exploring molecular landscapes using functionally annotated gene sets. however, pathway scores can also be used as predictive features in machine learning. that requires, firstly, increasing statistical power and biological relevance via a network enrichment analysis and, secondly, a fast and convenient procedure for rendering the original data into a space of pathway scores. however, previous implementations of nea involved multiple runs of network randomization and were therefore slow.  here, we present a new r package nearender which can transform raw'omics' features of experimental or clinical samples into matrices describing the same samples with many fewer nea-based pathway scores. this is done via a parametric estimation of the null binomial distribution and is thus much faster and less biased than randomization procedures. further, we compare estimates from these two alternative procedures and demonstrate that the summarization of individual genes to pathways increases the statistical power compared to both the default differential expression analysis on individual genes and the state-of-the-art gene set enrichment analysis. the package also contains functions for preparing input, modeling null distributions, and evaluating alternative versions of the global network.  beyond the state-of-the-art exploration of molecular data through pathway enrichment, score matrices produced by nearender can be used in larger bioinformatics pipelines as input for phenotype modeling, predicting disease outcomes etc. this approach is often more sensitive and robust than using the original data. the package nearender is complementary to the online nea tool evinet and, unlike of the latter, enables high performance of computations off-line. the r package nearender version is available at cran repository https://cran.r-project.org/web/packages/nearender/ electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords enrichmentnetwork analysisnetwork benchmarkr packagejoint 14th network tools and applications in biology workshop and 11th integrative bioinformatics symposium nettab/ib 2015bari, italy 14- october issue-copyright-statement© the author 2017 