BACKGROUND
gene regulation is one of the most important biological processes in living cells.
it is indispensable for adapting to changing environments, stimuli, and developmental stage and plays an essential role in the pathogenesis and course of diseases.
mechanistically, the transcription of dna into rna is predominantly controlled by a complex network of transcription factors  .
these proteins bind to enhancer or promoter regions adjacent to the genes they regulate  <cit> , which may enhance or inhibit the recruitment of rna polymerase and thereby activate or repress gene transcription  <cit> .
gene products also can be modified post-translationally via micrornas  degrading the transcript or inhibiting their translation  <cit> .
besides, a multitude of other mechanisms influence gene regulation, such as chromatin remodelling  <cit> , epigenetic effects  <cit> , and compound-building of transcription factors  <cit> .
distortion of regulatory processes is inflicted with various diseases  <cit> , especially with cancer  <cit> .fig.
<dig> transcription of dna into rna.
transcription factors  bind to distal or proximal tf binding sites  enhancing the binding of rna polymerase and activating the transcription of dna into rna
due to this importance, many efforts have been devoted to the elucidation of human regulatory relationships and networks.
wide-spread experimental techniques are transcriptome measurements to quantify gene and transcription factor co-expression  <cit> , chromatin immunoprecipitation  on chips or followed by sequencing for identifying binding patterns of specific tfs  <cit> , and bisulfite sequencing to find epigenetic signals of regulation  <cit> .
many large-scale datasets of such experiments have been published and are available in public repositories such as the gene expression omnibus   <cit> , the cancer genome atlas   <cit>  or the encyclopedia of dna elements   <cit> .
computational methods are also used, for instance, to identify transcription factor binding sites   <cit>  or to find known tfbs within the genome .
several databases have been created which store relevant information, such as lists of binding motifs  or targets of regulatory mirnas  <cit> .
such measurements and predictions are used by network reconstruction algorithms to predict regulatory relationships and regulatory networks  <cit> .
a plethora of different methods have been proposed, ranging from purely qualitative methods  <cit>  over simple statistical approaches  <cit>  to more advanced probabilistic frameworks  <cit> .
early methods were plagued by insufficient data and a general scarcity of background knowledge, which led to rather unstable results  <cit> .
this situation has changed dramatically over the last years, as results of more and more large screens have been made publicly available  <cit>  and also the knowledge on principal regulatory relationships has increased  <cit> .
this, in turn, has increased the interest in methods which predict genome-wide networks using a systematic, unified, mathematical framework.
here, we review five rather recent methods and conduct a quantitative comparison of their results with the goal to identify their mutual strengths and weaknesses.
they all have in common that they assume both the set of regulators  to be known and the topology of the regulatory network to be given.
by combining this background knowledge with specific omics data sets, especially transcriptome data, they try to infer the activity of regulators in a certain experimental condition or disease using mathematical optimization.
all presented methods are global methods in the sense that they compute activities genome-wide , thus removing the shortcomings of local methods which ignore cross-talk between sub-models and global effects within samples.
the methods predominantly produce a ranked list of regulators, sorted by their activity in a given group of samples; given that a multitude of biological influences is ignored during inference, especially kinetic and temporal effects, their goal cannot be to produce absolute snapshots of regulatory activity.
we describe each method in detail and compare them with respect to the most important properties, such as the data being used, the method applied for deriving optimized activity values, or the evaluation performed to show effectiveness.
we further implemented a quantitative comparison including four of the presented methods to objectively analyze their results.
as contrast, we also include aracne  <cit>  as sixth method; this algorithm uses only local reasoning and requires no background knowledge, but is still rather popular.
methods
we describe in detail five methods which infer transcription factor activity from omics data sets using a background network of transcription factors and the genes they regulate.
all use some form of mathematical optimization.
to emphasize the common ground of these at-first-sight rather different methods, we explain their underlying models using a simple framework for defining the relationships of transcription factors and genes.
this framework is presented first; it should be understood as a least common denominator, not as a proper method for network inference by itself.
we then describe five recently published methods for genome-wide tf activity estimation as extensions or constraints to this general framework, namely the approach by schacht et al.
<cit>  , racer  <cit> , rabit  <cit> , ismara  <cit>  and birte  <cit> .
additionally, we contrast these more comprehensive methods with the local inference algorithm aracne  <cit> , a popular tool for the de-novo reconstruction of gene regulatory networks.
key properties of all methods  are summarized in table  <dig> table  <dig> overview of methods for estimating regulatory activity from transcriptome data comparing input data, modelling, computational aspects and outcome variables
gi,s^=c+∑tβtbt,iθa,tactt,s+θg,tgt,s
with
actt,s=∑ibt,igi,s∑ibt,i,θa,t+θg,t= <dig> θa,t,θg,t∈01
- mixed-integer linear programming
t
- decision for each tf if θ
a,tor θ
- copy number variation
- dna methylation
- mirna expression signals
- tf binding information
1) gi,s^=c+θcnv,scnvi,s+θdm,sdmi,s+∑tβt,sbt,i+∑miβmi,sci,mimirnami,s
2) gi,s^=c˜+θ˜i,cnvcnvi,s+θ˜i,dmdmi,s+∑tγi,tβt,s+∑miγi,miβmi,s
t,s and β
mi,s
2) tf-gene γ
i,t and mirna-gene γ
- somatic mutations
- dna methylation
- copy number variation
- tf binding info
gt^=∑fθfbf,i+∑tβtbt,i
- annotation of promoters
- transcripts and associated promoters
gp,s^=cp+cs+∑mnp,mβm,s
- bayesian procedure, ridge regression
- gaussian prior for β
m,swith set of tfs and mirnas binding to sites of these motifs
- predicted target promoters, associated transcripts and genes
- network of known interactions between predicted targets and predicted regulatory interactions
- mirna, tf measurements, cnv
ld,θr=pd|r,θ=∏d^pd^|r,θ=∏d^∏c∏ipd^ic|rc,θ
gene expression data is named “g” with index i, estimated parameters with “β”, tf binding information with “b”, tfs with “t”, samples with “s”, mirnas with “mi” and model constants with “c”.
other variables are explained in the text
mathematical framework
to combine regulatory networks and quantitative omics data and to thereby deduce regulatory activity, all methods described here use a genome-wide mathematical model.
sample specific gene expression values g
i,s, derived from one biological condition, i.e., grouped into a single class, for in total g genes and s samples need to be provided as input.
the background regulatory network is represented as a directed graph where the nodes designate regulators and regulated entities  and directed edges indicate a regulatory relationship between the two connected nodes, for example the influence of a tf on the expression of a gene .fig.
<dig> general scheme of a tf – gene network where all t tfs are connected to each other and can regulate all of the g genes
we will use the variable t for regulators, i for regulated entities, and b
t,i for the strength of an edge from a tf/mirna t to a gene i representing, for instance, a binding affinity.
as abstract framework for explaining the different methods we propose a simple linear model predicting gene expression gi,s^ of gene i in sample s in terms of the activity of all t transcription factors β
t,s, which regulate i, and the binding affinities b
t,i.
in contrast to fig.
<dig>  where tfs can influence each other, this model ignores tf – tf relations and feedback loops: gi,s^=∑t=1tβt,sbt,i
given this model and a set of quantitative measurements of gene expression g
i,s, the goal of the mathematical optimization is to find parameters β such that the sum of squared errors of measured vs predicted gene expression over all genes and samples is minimized using a certain norm, for example the l <dig> norm: min∑i=1g∑s=1sgi,s−gi,s^ <dig>
estimation of tf activity by the effect on their target genes  <cit>
the idea of this method is to use the expression levels of tf’s target genes to infer their integrated effect .
the method uses expression data and database curated tf binding information as input whereby the tf – gene network is restricted to genes regulated by more than  <dig> tfs and tfs with at least  <dig> target genes.
the model is closely related to the abovementioned general framework, only adding a term for the sample specific effect of a tf.
specifically, the activity of a tf is modelled linearly by its cumulative effect on its target genes normalized by the sum of target genes or the tf’s gene expression level:fig.
<dig> flow chart of the approach by schacht et al.
the input data sets  are partly filtered and passed to a linear regression model  which calculates an activity value for each tf
gi,s^=c+∑tβtbt,iθa,tactt,s+θg,tgt,s where gi,s^ denotes the predicted gene expression of gene i in sample s, c is an additive offset, β
t describes the estimated activity of tf t and b
t,i refers to the underlying strength of the relation between tf t and gene i reflecting the binding affinity.
the estimated effect of a tf in a certain sample is calculated via the switch-like term in parentheses, where either the activity definition actt,s=∑ibt,igi,s∑ibt,i or the gene expression of the tf itself g
t,s is taken into account using the restrictions θ
a,t, θ
g,t ∈ { <dig>  1} and θ
a,t + θ
g,t =  <dig>  this switch term represents a meta-parameter to find the best model and has no biological interpretation.
the model outputs an activity value and the information which switch parameter is chosen for each tf of the reduced network.
during the optimization, the sum of error terms  is minimized which is achieved via mixed-integer linear programming using the gurobi  <dig>  optimizer.
<dig> the authors of this method state that the activity definition  was used in 95% of their test cases, but the switch-like combination of both terms yielded still better optimization results.
in the paper, the optimization task is greatly simplified as the model is computed for each gene separately and allows only a maximum number of  <dig> regulating tfs.
the tf – gene network indicating the strength of a relation between a tf and a gene is created for  <dig> tfs using knowledge from the commercial metacore™ database, <dig> chea  <cit>  and encode  <cit> .
due to the restriction of the network mentioned above, the actual model is then based on  <dig> tfs and  <dig> target genes only.
evaluation of the results was performed using expression data from  <dig> cell lines of the nci- <dig> panel  <cit>  and from melanoma cell lines   <cit> .
a sample based leave-one-out and 10-fold cross validation of predicted and measured gene expression yielded pearson correlation scores of about  <dig>  for both data sets.
a gene set enrichment analysis of the target genes for tfs modelled by the activity definition yielded  <dig> significantly enriched concepts including cell cycle, immune response and cell growth for the data from the nci- <dig> panel.
additionally, a t-test was computed between melanoma and other cell lines of the nci- <dig> panel to find differentially expressed genes of melanogenesis.
for the resulting genes, regulation models were built and used to predict gene expression in the melanoma cell line data set yielding good prediction performances.
racer  <cit>
racer  aims to integrate generic cell-line data with sample-specific measurements using a two-stage regression .
firstly, sample-specific regulatory activities for tfs and mirnas are calculated.
subsequently, general tf/mirna – gene interactions are derived.fig.
<dig> scheme of racer method.
the input data sets  are passed to a two-step linear regression model  which calculates sample specific activity values for each regulator and determines the most predominant regulators
compared to our general framework, racer includes additionally mirna binding information.
it assumes a linear combination, which is not further justified, of the regulatory effects of tfs and mirnas on mrna level.
racer can incorporate a variety of sample specific data including mrna and mirna expression values, cnv and dna methylation.
optimization is applied twice to reduce model complexity, where the method first infers sample-specific tf and mirna activities and uses these, in a second step, to compute general tf/mirna – gene interactions.
in the first regression step, mrna, mirna, cnv and dna methylation data are used to calculate the sample specific activities: gi,s^=c+θcnv,scnvi,s+θdm,sdmi,s+∑tβt,sbt,i+∑miβmi,sci,mimirnami,s where gi,s^ denotes the predicted gene expression of gene i in sample s, c is an intercept, β
t,s describes the estimated activity of tf t in sample s and b
t,i is the tf – gene binding score for tf t and gene i. the parameter β
mi,s stands for the estimated activity of mirna mi in sample s and is multiplied by c
i,mi, the number of conserved target sites on 3’utr of the target gene i for mirna mi, and by the expression level of mirna mi in sample s. θ
cnv,s  are the regression parameters for cnv signals cnv
i,s .
using β
t,s and β
mi,s from the first regression step, tf – gene and mirna – gene interactions across all samples are calculated in a second model: gi,s^=c˜+θ˜i,cnvcnvi,s+θ˜i,dmdmi,s+∑tγi,tβt,s+∑miγi,miβmi,s where the sums apply only to a number of selected tfs and mirnas with nonzero binding signals b
t,i >  <dig> and conserved target sites c
i,mi >  <dig>  the resulting parameters γ
i,t and γ
i,mi indicate the strength of a tf/mirna – gene relationship across all samples.
to obtain robust estimates, γ
i,mi is additionally weighted by the averaged activities of the mirna.
in each of the two regression steps, the optimization criterion is to minimize the sum of squared errors with l <dig> penalty on the linear coefficients to induce a sparse solution and to set irrelevant parameters to zero after the fitting.
this sparse lasso solution is obtained through elastic-net regularized generalized linear models.
a supplementary feature selection procedure comparing the full model to a restricted model leaving one tf or mirna out provides the most predominant tf/mirna regulators.
tf binding scores are collected from the generic cell line of erythroleukemia cells k <dig> from encode for  <dig> tfs and  <dig> genes.
further, the number of conserved target sites on 3’utr is taken from sequence-based information from targetscan for  <dig> mirnas and  <dig> genes.
the racer method is implemented in r and publicly available under http://www.cs.utoronto.ca/~yueli/racer.html.
the method was evaluated using expression data from an acute myeloid leukemia  data set from tcga with  <dig> samples  <cit>  via a sample based 10-fold cross validation on the prediction of gene expression.
to assess the quality of predictions, the spearman rank correlation was calculated resulting in a reassuring value of approximately  <dig> .
further, the full model was compared to models excluding one type of the input variables.
the full model performed best and a substantial reduction of spearman correlation was observed by omitting tf regulation  and dna methylation .
racer also performed with competitive accuracy in predicting known mirna – mrna and tf – gene relationships compared to other methods like genmir++  <cit>  or encode tf binding scores  <cit>  using e.g., validated interactions from the mirtarbase  <cit>  and knockdown studies.
the feature selection procedure revealed  <dig> predominant transcriptional regulators in the aml dataset.
using their associated targets, a functional enrichment analysis showed that dna repair and the tumor necrosis factor pathway were enriched.
when applying this panel to cluster patients at different cytogenetic risks, the clustering pattern of the regulatory activities was largely consistent with the risk groups.
further, a literature survey on aml showed that many tf regulators among the top predictions had a role in leukemogenesis.
rabit  <cit>
regression analysis with background integration  is a method for finding expression regulators in cancer by a large scale analysis across diverse cancer types.
it integrates tf binding information with tumor profiling data to search for tfs driving tumor-specific gene expression patterns .
it can be applied to predict cancer-associated rna-binding protein  recognition motifs which are key components in the determination of mirna function  <cit> .fig.
<dig> flow chart of rabit method.
the input data sets  are passed to a linear regression model  which calculates sample specific activity values for each regulator and determines general regulatory activities
in contrast to our general framework, rabit can, like racer, make use of cnv and dna methylation data additionally integrating promoter cpg content and promoter degree information  and takes rbp or tf binding information as regulatory input.
the computational model consists of three steps .
first, rabit tests in each tumor whether the target genes, identified by the beta method  <cit> , show differential expression compared to the normal controls including a control for background effects from cnvs, promoter dna methylation, promoter cpg content and promoter degree: gi^=∑fθfbf,i+∑tβtbt,i where gi^ represents the predicted differential gene expression between tumor and normal samples in gene i, b includes values of the f different background factors for gene i, b contains rbp or tf binding information and θ and β are the respective regression parameter vectors.
the regression coefficients β are estimated by minimizing the squared difference between measured and predicted gene expression.
the regulatory activity score for each tf/rbp is defined by a t-value  and its significance by the corresponding t-test.
if multiple profiles exist for the same tf from different conditions or cell lines, the profile with the highest absolute value of tf regulatory activity score is selected.
in a second step, a stepwise forward selection is applied to find a subset of tfs among those screened in step one optimizing the model error.
lastly, tfs with insignificant cross-tumor correlation are removed from the results.
computationally, the regression coefficients are calculated via the efficient frisch-waugh-lovell method.
tf binding information is taken from  <dig> tf chip-seq profiles from encode representing  <dig> tfs and  <dig> cell types.
additionally, recognition motifs for  <dig> rbps and their putative targets are collected by searching recognition motifs over the 3’utr regions  <cit> .
an implementation of the rabit method can be downloaded from http://rabit.dfci.harvard.edu/download.
rabit was applied to  <dig> tumor profiles of  <dig> cancer types from tcga using gene expression, somatic mutation, cnv and dna methylation data.
to systematically assess the results, the cancer relevance level of a tf was calculated as percentage of tumors with the tf target genes differentially regulated .
a comparison to cancer gene databases, i.e., the nci cancer gene index project  <cit> , the bushman laboratory cancer driver gene list  <cit> , the cosmic somatic mutation catalog  <cit>  and the ccgd mouse cancer driver genes  <cit> , showed a consistent picture.
further, rabit’s performance was compared to other regression models like lar or lasso where rabit had the best classification results when classifying all tfs into three categories by nci cancer index and achieved better cross-validation error and shorter running time.
the regulatory activity of rbps showed that some alternative splicing factors could affect tumor-specific gene expression by binding to target gene 3’utr regions.
ismara  <cit>
in contrast to the previous three methods and to our general framework which directly scores tfs or other regulators, ismara  infers the activity of regulatory motifs  and thereby indirectly deduces the effects of tfs and mirnas .
ismara is a web service where no parameter settings or specific processing of the input data, gene expression or chip-seq data are necessary.
it can also be used to calculate regulatory activity differences between samples and consider replicates or data from time series.fig.
<dig> ismara model scheme.
the input data sets  are passed to a linear regression model  which calculates motif activities and determines associated regulators
ismara takes sample specific measurements and information about regulatory motifs for tfs and mirnas into account.
based on the input of gene expression data or chromatin state measurements, the input signal is calculated for each promoter in each sample.
the input signal levels are modelled linearly in terms of the binding site predictions and unknown motif activities: gp,s^=cp+cs+∑mnp,mβm,s where gp,s^ refers to the input signal for a promoter p in sample s, c
p
and c
s are intercepts for each promoter and sample, n
p,m summarizes the tf/mirna binding site predictions  and β
m,s stands for the estimated motif activities.
like in the other presented methods, the optimization criterion is to minimize the sum of squared error terms between predicted and measured gene expression.
primarily, ismara provides the inferred motif activity profiles  sorted by significance and a set of tfs and mirnas that bind to these motifs representing the key regulators.
further, a list containing their predicted target promoters, associated transcripts and genes, a network of known interaction between these targets and a list of enriched gene ontology categories is displayed.
the web service ismara is available under http://ismara.unibas.ch.
ismara employs a bayesian procedure with a gaussian likelihood model and a gaussian prior distribution for β
m,sto avoid overfitting.
information about regulatory motifs is provided via the annotation of promoters based on deep sequencing data of transcription start sites.
to obtain a set of promoters and their associated transcripts, the 5’ ends of mrna mappings from ucsc genome database are clustered with the promoters.
tf binding site predictions in the proximal promoter region are collected using  <dig> position weight matrices representing  <dig> tfs from jaspar, transfac, motifs from the literature and their own analyses of chip-seq and chip-chip data.
additionally, mirna target sites for about  <dig> seed families are annotated in the 3’utrs of transcripts associated with each promoter.
for evaluation, ismara was applied to data from well-studied systems and results were compared to the literature.
inferred motif activities were highly reproducible and even more robust than the expression profiles from which motif activities were derived.
when comparing samples from  <dig> human cell types  from younger and older donors, ismara was able to identify a key regulator of aging-related changes in expression of lysosomal genes.
a joint analysis of the human gnf atlas of  <dig> tissues and cell lines  <cit>  and the nci- <dig> reference cancer cell lines  <cit>  revealed that many of the top dysregulated motifs were well-known in cancer biology like hif1a and has-mir- <dig> mirna.
they also suggested novel predictions for regulating tfs in innate immunity, mucociliary differentiation and cancer.
birte  <cit>
birte  takes a mathematically different approach compared to the abovementioned methods integrating tf/mirna target gene predictions with sample specific expression data into a joint probabilistic framework .
compared to our general scheme of a tf – gene network , birte takes the tf/mirna – gene network without the interactions between regulators to estimate regulatory activities and infers the network between regulators in a second step.fig.
<dig> scheme of birte method.
the input data sets  are passed to a likelihood model  which determines active regulators
birte takes as input differential gene expression data , an underlying regulatory network including tf/mirna – target gene binding information and optionally cnv data, mirna and tf expression measurements.
as opposed to our general framework, birte defines a likelihood model for the set of active tfs/mirnas  based on the entire gene expression data d and certain model parameters θ: ld,θr=pd|r,θ=∏d^pd^|r,θ=∏d^∏c∏ipd^ic|rc,θ
here d represents the set of all available experimental data including mrna, cnv, mirna and tf expression data and d
ic refers to its ith feature measured under experimental condition c. the condition specific hidden state variables r
c are estimated with help of the markov chain monte carlo  method where a regulator can switch from an active to an inactive state  or an inactive and an active regulator exchange their activity states .
thereby, the posterior probability for each regulator and condition to influence the expression of its target genes is estimated.
simultaneously, a variable selection procedure is applied to achieve sparsity of the model.
the optimization goal is not, as one would expect, to return the configuration with highest posterior probability among all sampled ones but to take marginal selection frequencies during sampling into account and filter those above a defined cutoff.
after the determination of active regulators, the associated transcriptional network containing tfs and mirnas is inferred from the observable differential expression of target genes and target gene predictions for individual regulators.
in practice, the stochastic sampling scheme based on mcmc allows swap operations only when regulators show a significant overlap of regulated targets.
the variable selection procedure is implemented via a spike and slab prior  <cit>  which can integrate prior knowledge about the activity of regulators.
to infer the associated transcriptional network, nested effects model   <cit>  structure learning is applied.
an input mirna – gene network is constructed based on mirmap  <cit>  for  <dig> mirnas.
the tf – target gene network with  <dig> tfs is compiled by computing tf binding affinities to promoter sequences according to the trap model  <cit>  using data from ensembl, transfac, jaspar and metacore™.
an implementation of birte is available for r on bioconductor under https://bioconductor.org/packages/release/bioc/html/birte.html.
several simulations were conducted to study model behavior.
on the basis of a human regulatory sub-network and accordingly simulated expression data of  <dig> target genes birte was compared to birta  <cit> , gemula  <cit>  and a hypergeometric test and further to other network reconstruction algorithms like aracne  <cit> , genie <dig>  <cit>  and genenet  <cit> .
birte performed best in regulator activity predictions including a favorable computation time and was robust against false positive and false negative target gene predictions.
additionally, birte was applied to an e.coli growth control and to a prostate cancer data set including  <dig> normal and  <dig> cancer samples from geo  with corresponding array data from  <dig> human mirnas  and the results showed a principal agreement with the biological literature.
aracne  <cit>
we compare aracne   <cit>  as an established, yet local, tool for the reconstruction of gene regulatory networks to the previous five recent genome-wide approaches.
the algorithm is background knowledge-free and identifies transcriptional interactions based on mutual information including non-linear and non-monotonic relationships and distinguishes between direct and indirect relationships .
aracne is a free tool available under http://califano.c2b <dig> columbia.edu/aracne.fig.
<dig> aracne flow chart.
the input data set  is used to calculate pairwise mutual information where indirect interactions are removed  and which allow a reconstruction of the gene regulatory network
aracne uses as input only microarray expression profiles and estimates candidate interactions by calculating the pairwise gene expression profile mutual information i defined as i=ii,j=s+s−s where s denotes the entropy.
i
i,j measures the relatedness of genes g
i and g
j and equals zero if both are independent.
in a second step, the mutual information values are filtered using an appropriate threshold depending on the distribution of all mutual information values between random permutations of the original data set and indirect interactions are removed.
computationally, a gaussian kernel operator is used to calculate mutual information scores.
in a subsequent step, the data processing inequality   <cit>  is applied to remove probably indirect candidate interactions.
the dpi states that if the genes g
i and g
k interact only through a third gene g
j, then i≤min,i)
thus, the least of the three mutual information scores can come from indirect interactions only  <cit> .
aracne’s performance was evaluated on the reconstruction of realistic synthetic datasets  <cit>  and on an expression profile dataset consisting of about  <dig> b lymphocytes derived from normal, tumor-related and experimentally manipulated populations  <cit>  against relevance networks and bayesian networks.
regarding the synthetic networks, aracne had consistently better precision and recall values compared to the two other algorithms and reached very good precision at significant recall levels.
it recovers far more true connections and fewer false connections than the other methods with better performance on tree-like topologies compared to scale-free topologies.
a reconstructed b-cell specific regulatory network was found to be highly enriched in known c-myc targets where about 50% of the predicted genes to be first neighbors were reported in the literature.
RESULTS
we described five recent methods for the genome-wide inference of regulatory activity, namely the approach by schacht et al., racer, rabit, ismara, and birte.
they all assume the topology of the regulatory network to be known, cast activity estimation as an optimization problem regarding the difference between predicted and measured values, take different types of sample specific omics data into account, and eventually produce a list of regulators like transcription factors or mirnas, ranked by their estimated activities in the samples under study.
we also included aracne which is background knowledge-free and uses only local dependency measures to reconstruct a regulatory network and indirectly infer activities.
all of the presented methods essentially follow the same goal, i.e., accurate ranking of regulatory activity, but differ in the types of measurements being integrated, the background knowledge necessary for their application, the complexity and refinement of the underlying model of gene regulation, and the concrete paradigm used for solving the optimization problem.
most of the methods, except for the approach by schacht et al., are available online via a downloadable implementation, a web service, or an r package providing an operable solution for the interested user.
whereas an overview of the main features of each method ca be found in table  <dig>  we now first compare the algorithms regarding their general properties in a descriptive way.
the data sets used for evaluation vary between all methods.
therefore, we further implemented an evaluation framework to compare the method by schacht et al., racer, rabit and birte in an objective and quantitative way.
we used experimental data of three publicly available data sets from tcga  <cit>  and a regulatory network as background knowledge.
we first used only mrna expression data as input to the four methods to ensure the result’s comparability, whereas in a second evaluation step, also other omics data sets were included where possible.
we further analyzed the relevance of regulators found by different methods using a literature search.
general properties
experimental data types included
the methods differ in the types of measurements being integrated, which corresponds to the level of detail of their model of gene regulations.
all six methods use mrna as input.
racer, rabit and birte can also integrate cnv, dna methylation, tf/mirna expression data, or somatic mutations.
ismara calculates an input signal from microarray, rna-seq, or chip-seq data.
additionally, all presented methods use prior knowledge about the underlying regulatory network.
these networks are extracted from different data sources and pre-processed in different manners.
all methods require at least knowledge about tf – gene relationships, yet racer, birte and ismara also incorporate information about mirnas.
when using rabit, the user can choose whether to provide tf or rna binding protein information.
the approach of schacht et al.
and birte extract regulatory information partly from the commercial metacore™ database, whereas the other methods use only publicly available databases, like encode, jaspar or transfac.
the networks which are used for the evaluations published in the respective papers are publicly available for the case of racer , rabit  and birte .
neither schacht et al.
nor ismara make this data available.
mathematical models of regulatory activity
the methods use different mathematical models to infer regulatory activity.
the approach by schacht et al., racer, rabit and ismara use linear regression whereas birte applies a probabilistic framework.
aracne, as a local method, is based on mutual information.
racer and rabit can be seen as extensions of the approach by schacht et al.
since they essentially use the same model structure but incorporate more input data types and more classes of regulatory information.
further, racer applies a two-stage regression to infer regulatory activity.
optimization frameworks
for assessing regulator activities, schacht et al., racer, rabit and ismara minimize the sum of error terms between measured and predicted gene expression.
however, the methods use rather different algorithms for solving the resulting optimization problem, and also apply different constraints to achieve model sparsity, robustness of inference, and feature selection.
in the approach by schacht et al., the regression model is computed for each gene separately and allows only a maximum number of six regulating tfs.
racer uses a lasso approach, while ismara follows a bayesian model that infers regulator activities as posterior distributions.
lasso can be interpreted as a bayesian model using laplacian priors instead of gaussian priors in the regression framework obtaining point estimates of the regulatory activities and enforcing sparseness of the solution  <cit> .
in contrast, birte uses a likelihood model with a spike and slab prior to induce model sparsity.
this approach implements a selective shrinkage of model coefficients such that estimates are less biased compared to a lasso prior  <cit> .
with the help of the spike and slab prior, sparsity can be controlled in a variable dependent manner allowing the inclusion of prior belief in the activity of each regulator  <cit> .
computed outputs
schacht et al.
and birte determine activity of regulators over all samples at once, whereas racer and birte first infer sample-specific activities which are combined to cross-tumor activities only in a second optimization step.
in contrast, ismara in first place infers motifs activity; these activities are used to deduce the effects of tfs and mirnas by their motif binding profiles.
ismara primarily provides sample specific tf and mirna activity but also offers an option to group samples and compare average regulatory activity between different conditions.
like birte and aracne, it also infers the network of the regulators themselves.
methods and data sets used for evaluation
the type and extent of evaluation performed for the different methods vary greatly.
they range from direct application to biological problems over the comparison of results to the biological literature to simulation studies.
all methods published evaluations results on publicly available datasets, e.g., from the national cancer institute, tcga or geo, but unfortunately address different tissues and cancer types.
sample-based cross-validation is applied in the work by schacht et al., racer, rabit and ismara.
the first two of these methods use correlation coefficients between measured and predicted gene expression for assessing prediction quality.
racer, rabit and birte compare their results to the outcome of other algorithms and to those of restricted models, for example excluding one type of the input variables.
all methods search the literature to compare their predictions to previously published studies on the respective biological question.
overall, ismara provide the most extensive biological evaluation using a battery of relevant use cases, whereas birte excels in systematic simulation studies.
sadly, there are very few works which compare any of the methods presented on the same problem; the only result we are aware of compared aracne and birte regarding their performance in network reconstruction on simulated data, in which birte attained higher robustness against false positive and false negative target gene predictions  <cit> .
quantitative comparison
although certain evaluation steps were carried out for all methods, results in the original papers are not comparable as they used different input datasets, different background regulatory networks, and different evaluation metrics.
therefore, in addition to the comparison of general properties of the methods, we implemented an evaluation framework using three independent and publicly available test data sets to compare the method by schacht et al., racer, rabit and birte in an objective and quantitative way.
all evaluated methods were given the same regulatory network as input.
data sets
for the evaluation we used experimental data from tcga  <cit>  for three cancer types: colon adenocarcinoma , liver hepatocellular carcinoma  and pancreatic adenocarcinoma .
for all three cancer types, mrna expression, cnv, dna methylation and mirna expression data is available for primary tumor and normal tissue samples.
these data sets are openly accessible via the nci genomic data commons data portal <dig> or the nci genomic data commons legacy archive <dig> .
for mrna gene expression we used processed rna-seq data in the form of fpkm  values.
the files included ensembl gene ids which were converted to hgnc symbols using the ensembl  <cit>  biomart tool <dig> to match the ids of the tf – gene network.
in two cases, when multiple ensembl gene ids mapped to one hgnc symbol, we chose the gene with highest log <dig> fold change between case and control group.
mirna expression was given as rpm  measurements.
both mrna and mirna data were centered using a weighted mean such that the mean of the case group equaled the negative mean of the control group, and normalized via a weighted standard deviation.
cnv data was retrieved as masked copy number segment where the y chromosome and probe sets with frequent germline copy-number variation had already been removed.
chromosomal regions were mapped to genes using the r package biomart  <cit> .
if multiple records mapped to one gene, the median of the segment mean values was calculated.
for dna methylation data we used the beta-values of illumina human methylation  <dig> arrays as methylation scores.
multiple scores for the same gene were averaged within a sample.
we restricted our analyses to the samples for which all four input data types were available.
when multiple measurements for one sample and data type were available, we used only the first one in alphabetical order of the file name.
after this selection procedure,  <dig> samples remained for coad,  <dig> for lihc and  <dig> for paad.
a list including sample and file information is available in additional file  <dig>
together with the experimental data, all evaluated methods were given the same regulatory network as input.
we used a publicly available human tf – gene network  <cit>  based on a text-mining approach and complemented it with tf – gene interactions from the public transfac <dig> database  <cit> .
this network included  <dig> interactions between  <dig> tfs and  <dig> genes.
the network is provided in additional file  <dig>
evaluated methods
we conducted the quantitative comparison for the method proposed by schacht et al., racer, rabit and birte.
ismara was not included since it is  only available as a web service,  can only be used with its own, proprietary underlying regulatory network model, and  requires the upload of raw data which is prohibited by tcga’s terms of use.
also aracne  <cit>  was not included in the quantitative evaluation since it does not use background knowledge and we therefore consider its results as incomparable to the other methods.for the approach by schacht et al.
we re-implemented their method as closely as possible to the original design using python and the cuneiform workflow language  <cit> .
due to the high number of integer parameters in the original method, the complexity of optimizing the whole network at once would have by far exceeded computational measures.
therefore, like in the original paper, we computed the model for each gene separately and restricted the number of regulating tfs per gene to six.
we added a second step where we used these tf – gene interactions building a sub-network to optimize tf activity globally to describe the interplay of the tfs’ effects on their target genes.
as in the implementation of schacht et al., we used the gurobi optimizer.7
for racer we used the available r scripts <dig> and extracted the resulting sample-specific regulatory activities.
rabit published a c++ implementation which they provide on their website <dig> and which we used with the fdr option set to  <dig>  as rabit takes differential expression into account, we used the difference of expression values between case and control group as input and ordered the tfs by t-value as proposed in the rabit paper.
birte is available as a bioconductor r package.
we used r version  <dig> .
<dig> with birte version  <dig> .
<dig> and applied the method “birtelimma” to estimate regulatory activities with the options niter and nburnin set to  <dig>  as birte has a randomized component, the resulting tf activities are not exactly the same for different runs.
we averaged the final activity scores over  <dig> iterations of birtelimma.
for our re-implemented method by schacht et al.
and racer we computed separate models for case and control group and ranked the tfs by their activity difference between the two groups.
to ensure the result’s comparability, we first used only mrna expression data as input to the four methods.
in a second evaluation, we included also other omics data sets where possible.
birte was evaluated on mrna and cnv data, rabit on mrna, cnv and dna methylation data, and racer additionally used mirna expression as input.
we obtained lists with the regulators ranked according to the absolute value of their computed activity for each cancer type and method, with and without the use of additional inputs.
for each cancer type we calculated the size of the overlaps in the four different results using the top  <dig> and top  <dig> regulators.
the results for the top  <dig> regulators using either only mrna or multiple omics data sets as input are shown in table  <dig> table  <dig> hgnc symbols of the top  <dig> regulators found by each method for coad , lihc  and paad  and the use of only mrna data as input  and multiple input data sets .
tfs with equal activity values are marked with*.
tfs found by several method’s top  <dig> are marked in bold , blue , red  or yellow
only mrna as input
when only mrna is used as input, one tf is commonly found by the three methods racer, rabit and birte in each data set, respectively: phox2b for coad, epas <dig> for lihc and elf <dig> for paad.
a literature search of these tfs and their targets revealed clear associations to the respective cancer type.
the tf obtained commonly for coad, phox2b, is related to tlx <dig>  a gene which has been shown to play a role in the tumorigenesis of gastrointestinal stromal tumors  <cit> .
epas <dig>  which was found in the lihc top  <dig> tfs of three methods, is linked to cxcl <dig>  which plays an important role in metastasis formation of hepatocellular carcinoma by promoting the migration of tumor cells  <cit> .
for paad, three methods ranked tf elf <dig> high, which is related to  <dig> genes in our network, inter alia to brca <dig> and lyn.
mutations in the brca <dig> gene have been implicated in pancreatic cancer susceptibility  <cit> , whereas the knockdown of lyn reduced human pancreatic cancer cell proliferation, migration, and invasion  <cit> .
these results underline that the methods are able to find biologically relevant information about regulation processes in cancer.
several tfs in the top  <dig> are found by two of the four methods for instance, racer and rabit have four common top  <dig> tfs  in the coad data set.
however, the top  <dig> tfs found by the method by schacht et al.
do not overlap with any top  <dig> tfs of the other methods in any data set.
the agreement of racer, rabit and birte in the top  <dig> tfs hints to the biological importance of the found tfs since this overlap is statistical significant as the probability of finding common tfs in three sets of ten randomly chosen ones out of  <dig> tfs  is below  <dig> .
additionally, the methods do identify different tfs for different data sets, indicating the importance of the actual cancer specific mrna expression values and that results are not dictated by the background network.
the results for the number of overlapping regulators in the top  <dig> between the four methods and the three different data sets are shown in fig.
<dig>  for rabit, only  <dig> tfs for coad  could be ranked since all other tfs had an activity value equal to zero.fig.
<dig> number of overlapping tfs in the top  <dig> of ranked tfs per method  in coad/lihc/paad is shown)
when looking at the overlap of three of the four methods, the number of overlapping tfs is still the highest for the triplet racer, rabit and birte.
for the lihc dataset two tfs are found in the top  <dig> of all four methods .
e2f <dig> is a downstream target of zbtb <dig>  which was associated to the expression of cell cycle-associated genes in liver cancer cells  <cit> .
two target genes of e2f <dig>  cdk <dig> and tp <dig> were also involved in liver cancer development  <cit>  and proposed as prognostic marker of poor patient survival prognosis in hepatocellular carcinoma  <cit> .
further, epigenetic alterations of the ednrb gene, a target of sox <dig>  might play an important role in the pathogenesis of hepatocellular carcinoma  <cit> .
even if the result of four methods finding two common tfs is not statistically significant , their association to liver hepatocellular carcinoma shows that the methods reach their goal of identifying relevant tfs.
however, when comparing different data sets, the methods tend to rank the same tfs under the top  <dig> to a greater or lesser extent.
for example, the overlap of all top  <dig> tfs of the three cancer types is only one tf for rabit and nine tfs for birte, but  <dig> tfs for the method by schacht et al.
and even  <dig> tfs for racer.
therefore, the results from rabit and birte seem to be more cancer type specific and less dependent on the regulatory network than the results from racer.
however, we did not specifically investigate the influence of the underlying network and its topology on the results which would be an interesting point for further research.
multi-omics data as input
when not only taking mrna into account but also mirna, cnv and dna methylation, the results are more difficult to compare between the methods, since they all use a different way of combining different types of data due to their models and implementations.
we are aware of the lower level of comparability of this approach regarding the multi-omics results in contrast to a scenario, where all methods are evaluated on the same set of input data.
however, we intended to use maximum set of input data for each method to cover the effect of the use of multiple omics data sets compared to only mrna as input.
birte was evaluated on mrna and cnv data, rabit on mrna, cnv and dna methylation data, and racer additionally used mirna expression as input.
whereas racer and rabit considered cnv or dna methylation data as one background factor and compute only one activity value, birte evaluated the influence of each cnv separately.
the results  show that racer exclusively ranks mirnas high; not a single tf is found among the top  <dig> regulators.
also, the influence of cnvs was high in lihc and paad.
however, the tfs that racer found in the top  <dig> when using only mrna data as input are still ranked high in the multi-omics scenario, e. g the coad top three tfs of the mrna results are ranked 13th, 16th and 14th in the results of the multi-omics input.
the difference of the results coming from the two input types is less for rabit: seven tfs are still in the top  <dig> for coad  when using cnv and dna methylation additionally to mrna data.
therefore, the contribution of additional input data seems not to be crucial for the performance of rabit.
birte considers each cnv as a potential regulator which increases the total number of regulators enormously.
still, two commonly present tfs in the top  <dig> of the coad data set  are found by either the sole mrna input and the multi-omics approach.
the overlap of the top  <dig> of rabit and birte in the multi omics case is considerable with three tfs in lihc , and one tf in paad .
three of them  were already found when using only mrna data as input.
the results for the use of different input data sets show that the top ranked regulators are drastically changed when using additionally mirna data in racer, but change less when only cnv or dna methylation data is provided in rabit and birte.
however, the results from multi omics analyses are difficult to compare since the combination of input data sets is not consistent across the three different methods.
discussion
background networks
a crucial input to the models is the underlying regulatory network which is needed to reduce the search space for actual regulatory activity.
however, the construction of comprehensive tf/mirna – gene regulatory networks is difficult for various reasons.
firstly, a comprehensive characterization of the human regulatory repertoire is lacking since only about half of the estimated  <dig> – <dig>  tfs in the mammalian genome is known  <cit> .
chip experiments, prone to a high false positive rate  <cit> , were used to identify tf binding patterns but each assay is limited to the detection of one tf in one condition and therefore tf binding has not been characterized for many tfs in most cell types.
further, the local proximity of a binding site to the transcriptional start site of a gene does not automatically implicate transcriptional regulation.
with regard to posttranscriptional regulators, the functions for only a few of the around  <dig>  different mirnas have been experimentally determined and current data on mirna targets is mostly based on computational predictions  <cit> .
generally, the knowledge about tf and mirna binding is scattered over the biological literature and different, partly commercial, databases, impeding the construction of comprehensive networks  <cit> .
therefore, any comparative evaluation of the methods presented here would have to make sure that the same background network is used for each computation.
besides, studies on the impact of network incompleteness or different error rates in networks would be important to assess the ability of the methods to cope with such common problems.
simulation studies will be vital in this regard.
the graph view on regulation
the modelling of regulatory networks as graphs, as used in all presented methods, is perhaps not the optimal representation for the underlying biological regulatory processes.
a graph cannot easily account for important effects such as tf complex formation and temporal and spatial synchronization of activities.
furthermore, tf binding is affected by chromatin state and the impact of posttranslational modifications on transcriptional activity which are difficult to include in a graph view on regulation.
the model’s dependence on the topological structure and the robustness to changes in the underlying network have not been evaluated or discussed in any of the presented methods even if these issues are known to have a severe influence in network analysis  <cit> .
underlying mathematical model
linear models, widely spread in different fields of science, provide a simple and easily understandable design but over-simplify the underlying biological processes.
nonlinear behavior, e.g., saturation effects, cannot be represented.
considering that the number of available samples is typically relatively small, the incorporation of many different data types and according parameters into the model could result in excessively complex designs prone to overfitting, but this issue lacks general awareness.
only two of the presented methods incorporate parameter priors , and two apply cross validation techniques to estimate prediction performance .
further, the effect of temporal buffering between tf binding and the actual effect on gene expression is not included in any of the methods.
comparability
all methods produce a ranked list of regulators.
comparing these results across different methods, even when applied on the same data set and using the same background network, is difficult since no generally accepted benchmarks are available.
therefore, there currently is no objective measure to designate a best method.
the closest comparable evaluation effort we are aware of is implemented in the “dream <dig> – network inference” challenge  <cit> , which targets gene regulatory network reconstruction.
the invited participants reverse-engineered a network from gene expression data, including a simulated network, and evaluated the results on a subset of known interactions or the known network for the in-silico case.
the approach of genie <dig>  <cit>  which trains a random forest to predict target gene expression performed best and the integration of predictions from multiple inference methods showed robust and high performance across diverse data sets.
however, an extensive competitive evaluation to determine active regulators based on a given regulatory network has, to the best of our knowledge not been carried out yet.
we therefore compared the results of four methods in a quantitative way.
the experimental data and the regulatory network we used as input are publicly available to ensure transparency of our results.
the results suggest that the methods are able to find biologically relevant information about regulation processes in cancer.
however, the result overlaps are rather low .
this seems surprising as all methods essentially follow the same goal, i.e., identification of the most differentially active tfs or genes.
we think further research is necessary to exactly characterize the specific strengths of each method.
furthermore, we did not investigate the influence of the underlying network on the results, which is another topic for further research.
CONCLUSIONS
despite their often rather involved procedures and models, none of the presented methods adequately reflects the biological reality of regulatory activity in cells.
a specific disease phenotype is rarely caused by a single gene but rather a product of the interplay of genetic variability, epigenetic modifications and post-transcriptional regulation mechanisms  <cit> .
the presented methods ignore a multitude of such factors like the effects of chromatin state and alternative splicing, nonlinear relationships between regulatory activity and gene expression, or kinetic and temporal effects.
furthermore, tfs themselves regulate the expression of other tfs forming feedback loops which are not considered in any of the presented methods.
nevertheless, the methods apparently are able to detect strong signals and produced promising results in terms of ranking transcription factors by their activity and are thus valuable tools for identifying biomarkers for specific phenotypes.
additional files
additional file 1: lists information about the samples and files from tcga included in our quantitative evaluation for all three cancer types .
additional file 2: includes an adjacency list of the connected nodes of the tf – gene network.
the list includes three columns  where each row indicates an association with the value of “edge” between a tf and a gene.
complexes of tfs are indicated with a separating “.” between their components.
<dig>
http://www.gurobi.com/products/gurobi-optimizer
<dig>
http://lsresearch.thomsonreuters.com/pages/solutions/1/metacore
<dig>
https://gdc-portal.nci.nih.gov
<dig>
https://gdc-portal.nci.nih.gov/legacy-archive
<dig>
http://www.ensembl.org/biomart/martview, release 87
<dig>
http://www.gene-regulation.com/pub/databases.html, release  <dig>
<dig> version  <dig> , available under a free academic license
<dig>
http://www.cs.utoronto.ca/~yueli/racer.html
<dig>
http://rabit.dfci.harvard.edu

 gene regulation is one of the most important biological processes in living cells. it is indispensable for adapting to changing environments, stimuli, and developmental stage and plays an essential role in the pathogenesis and course of diseases. mechanistically, the transcription of dna into rna is predominantly controlled by a complex network of transcription factors. these proteins bind to enhancer or promoter regions adjacent to the genes they regulate, which may enhance or inhibit the recruitment of rna polymerase and thereby activate or repress gene transcription. gene products also can be modified post-translationally via micrornas degrading the transcript or inhibiting their translation. besides, a multitude of other mechanisms influence gene regulation, such as chromatin remodelling, epigenetic effects, and compound-building of transcription factors. distortion of regulatory processes is inflicted with various diseases, especially with cancer.fig. transcription of dna into rna. transcription factors bind to distal or proximal tf binding sites enhancing the binding of rna polymerase and activating the transcription of dna into rna due to this importance, many efforts have been devoted to the elucidation of human regulatory relationships and networks. wide-spread experimental techniques are transcriptome measurements to quantify gene and transcription factor co-expression, chromatin immunoprecipitation on chips or followed by sequencing for identifying binding patterns of specific tfs, and bisulfite sequencing to find epigenetic signals of regulation. many large-scale datasets of such experiments have been published and are available in public repositories such as the gene expression omnibus, the cancer genome atlas or the encyclopedia of dna elements. computational methods are also used, for instance, to identify transcription factor binding sites or to find known tfbs within the genome. several databases have been created which store relevant information, such as lists of binding motifs or targets of regulatory mirnas. such measurements and predictions are used by network reconstruction algorithms to predict regulatory relationships and regulatory networks. a plethora of different methods have been proposed, ranging from purely qualitative methods over simple statistical approaches to more advanced probabilistic frameworks. early methods were plagued by insufficient data and a general scarcity of background knowledge, which led to rather unstable results. this situation has changed dramatically over the last years, as results of more and more large screens have been made publicly available and also the knowledge on principal regulatory relationships has increased. this, in turn, has increased the interest in methods which predict genome-wide networks using a systematic, unified, mathematical framework. here, we review five rather recent methods and conduct a quantitative comparison of their results with the goal to identify their mutual strengths and weaknesses. they all have in common that they assume both the set of regulators to be known and the topology of the regulatory network to be given. by combining this background knowledge with specific omics data sets, especially transcriptome data, they try to infer the activity of regulators in a certain experimental condition or disease using mathematical optimization. all presented methods are global methods in the sense that they compute activities genome-wide, thus removing the shortcomings of local methods which ignore cross-talk between sub-models and global effects within samples. the methods predominantly produce a ranked list of regulators, sorted by their activity in a given group of samples; given that a multitude of biological influences is ignored during inference, especially kinetic and temporal effects, their goal cannot be to produce absolute snapshots of regulatory activity. we describe each method in detail and compare them with respect to the most important properties, such as the data being used, the method applied for deriving optimized activity values, or the evaluation performed to show effectiveness. we further implemented a quantitative comparison including four of the presented methods to objectively analyze their results. as contrast, we also include aracne as sixth method; this algorithm uses only local reasoning and requires no background knowledge, but is still rather popular. methods we describe in detail five methods which infer transcription factor activity from omics data sets using a background network of transcription factors and the genes they regulate. all use some form of mathematical optimization. to emphasize the common ground of these at-first-sight rather different methods, we explain their underlying models using a simple framework for defining the relationships of transcription factors and genes. this framework is presented first; it should be understood as a least common denominator, not as a proper method for network inference by itself. we then describe five recently published methods for genome-wide tf activity estimation as extensions or constraints to this general framework, namely the approach by schacht et al., racer, rabit, ismara and birte. additionally, we contrast these more comprehensive methods with the local inference algorithm aracne, a popular tool for the de-novo reconstruction of gene regulatory networks. key properties of all methods are summarized in table  table overview of methods for estimating regulatory activity from transcriptome data comparing input data, modelling, computational aspects and outcome variables gi,s^=c+∑tβtbt,iθa,tactt,s+θg,tgt,s with actt,s=∑ibt,igi,s∑ibt,i,θa,t+θg,t= θa,t,θg,t∈01 - mixed-integer linear programming t - decision for each tf if θ a,tor θ - copy number variation - dna methylation - mirna expression signals - tf binding information 1) gi,s^=c+θcnv,scnvi,s+θdm,sdmi,s+∑tβt,sbt,i+∑miβmi,sci,mimirnami,s 2) gi,s^=c˜+θ˜i,cnvcnvi,s+θ˜i,dmdmi,s+∑tγi,tβt,s+∑miγi,miβmi,s t,s and β mi,s 2) tf-gene γ i,t and mirna-gene γ - somatic mutations - dna methylation - copy number variation - tf binding info gt^=∑fθfbf,i+∑tβtbt,i - annotation of promoters - transcripts and associated promoters gp,s^=cp+cs+∑mnp,mβm,s - bayesian procedure, ridge regression - gaussian prior for β m,swith set of tfs and mirnas binding to sites of these motifs - predicted target promoters, associated transcripts and genes - network of known interactions between predicted targets and predicted regulatory interactions - mirna, tf measurements, cnv ld,θr=pd|r,θ=∏d^pd^|r,θ=∏d^∏c∏ipd^ic|rc,θ gene expression data is named “g” with index i, estimated parameters with “β”, tf binding information with “b”, tfs with “t”, samples with “s”, mirnas with “mi” and model constants with “c”. other variables are explained in the text mathematical framework to combine regulatory networks and quantitative omics data and to thereby deduce regulatory activity, all methods described here use a genome-wide mathematical model. sample specific gene expression values g i,s, derived from one biological condition, i.e., grouped into a single class, for in total g genes and s samples need to be provided as input. the background regulatory network is represented as a directed graph where the nodes designate regulators and regulated entities and directed edges indicate a regulatory relationship between the two connected nodes, for example the influence of a tf on the expression of a gene.fig. general scheme of a tf – gene network where all t tfs are connected to each other and can regulate all of the g genes we will use the variable t for regulators, i for regulated entities, and b t,i for the strength of an edge from a tf/mirna t to a gene i representing, for instance, a binding affinity. as abstract framework for explaining the different methods we propose a simple linear model predicting gene expression gi,s^ of gene i in sample s in terms of the activity of all t transcription factors β t,s, which regulate i, and the binding affinities b t,i. in contrast to fig.  where tfs can influence each other, this model ignores tf – tf relations and feedback loops: gi,s^=∑t=1tβt,sbt,i given this model and a set of quantitative measurements of gene expression g i,s, the goal of the mathematical optimization is to find parameters β such that the sum of squared errors of measured vs predicted gene expression over all genes and samples is minimized using a certain norm, for example the l norm: min∑i=1g∑s=1sgi,s−gi,s^ estimation of tf activity by the effect on their target genes the idea of this method is to use the expression levels of tf’s target genes to infer their integrated effect. the method uses expression data and database curated tf binding information as input whereby the tf – gene network is restricted to genes regulated by more than tfs and tfs with at least target genes. the model is closely related to the abovementioned general framework, only adding a term for the sample specific effect of a tf. specifically, the activity of a tf is modelled linearly by its cumulative effect on its target genes normalized by the sum of target genes or the tf’s gene expression level:fig. flow chart of the approach by schacht et al. the input data sets are partly filtered and passed to a linear regression model which calculates an activity value for each tf gi,s^=c+∑tβtbt,iθa,tactt,s+θg,tgt,s where gi,s^ denotes the predicted gene expression of gene i in sample s, c is an additive offset, β t describes the estimated activity of tf t and b t,i refers to the underlying strength of the relation between tf t and gene i reflecting the binding affinity. the estimated effect of a tf in a certain sample is calculated via the switch-like term in parentheses, where either the activity definition actt,s=∑ibt,igi,s∑ibt,i or the gene expression of the tf itself g t,s is taken into account using the restrictions θ a,t, θ g,t ∈ {  1} and θ a,t + θ g,t =  this switch term represents a meta-parameter to find the best model and has no biological interpretation. the model outputs an activity value and the information which switch parameter is chosen for each tf of the reduced network. during the optimization, the sum of error terms is minimized which is achieved via mixed-integer linear programming using the gurobi optimizer. the authors of this method state that the activity definition was used in 95% of their test cases, but the switch-like combination of both terms yielded still better optimization results. in the paper, the optimization task is greatly simplified as the model is computed for each gene separately and allows only a maximum number of regulating tfs. the tf – gene network indicating the strength of a relation between a tf and a gene is created for tfs using knowledge from the commercial metacore™ database, chea and encode. due to the restriction of the network mentioned above, the actual model is then based on tfs and target genes only. evaluation of the results was performed using expression data from cell lines of the nci- panel and from melanoma cell lines. a sample based leave-one-out and 10-fold cross validation of predicted and measured gene expression yielded pearson correlation scores of about for both data sets. a gene set enrichment analysis of the target genes for tfs modelled by the activity definition yielded significantly enriched concepts including cell cycle, immune response and cell growth for the data from the nci- panel. additionally, a t-test was computed between melanoma and other cell lines of the nci- panel to find differentially expressed genes of melanogenesis. for the resulting genes, regulation models were built and used to predict gene expression in the melanoma cell line data set yielding good prediction performances. racer racer aims to integrate generic cell-line data with sample-specific measurements using a two-stage regression. firstly, sample-specific regulatory activities for tfs and mirnas are calculated. subsequently, general tf/mirna – gene interactions are derived.fig. scheme of racer method. the input data sets are passed to a two-step linear regression model which calculates sample specific activity values for each regulator and determines the most predominant regulators compared to our general framework, racer includes additionally mirna binding information. it assumes a linear combination, which is not further justified, of the regulatory effects of tfs and mirnas on mrna level. racer can incorporate a variety of sample specific data including mrna and mirna expression values, cnv and dna methylation. optimization is applied twice to reduce model complexity, where the method first infers sample-specific tf and mirna activities and uses these, in a second step, to compute general tf/mirna – gene interactions. in the first regression step, mrna, mirna, cnv and dna methylation data are used to calculate the sample specific activities: gi,s^=c+θcnv,scnvi,s+θdm,sdmi,s+∑tβt,sbt,i+∑miβmi,sci,mimirnami,s where gi,s^ denotes the predicted gene expression of gene i in sample s, c is an intercept, β t,s describes the estimated activity of tf t in sample s and b t,i is the tf – gene binding score for tf t and gene i. the parameter β mi,s stands for the estimated activity of mirna mi in sample s and is multiplied by c i,mi, the number of conserved target sites on 3’utr of the target gene i for mirna mi, and by the expression level of mirna mi in sample s. θ cnv,s are the regression parameters for cnv signals cnv i,s. using β t,s and β mi,s from the first regression step, tf – gene and mirna – gene interactions across all samples are calculated in a second model: gi,s^=c˜+θ˜i,cnvcnvi,s+θ˜i,dmdmi,s+∑tγi,tβt,s+∑miγi,miβmi,s where the sums apply only to a number of selected tfs and mirnas with nonzero binding signals b t,i >  and conserved target sites c i,mi >  the resulting parameters γ i,t and γ i,mi indicate the strength of a tf/mirna – gene relationship across all samples. to obtain robust estimates, γ i,mi is additionally weighted by the averaged activities of the mirna. in each of the two regression steps, the optimization criterion is to minimize the sum of squared errors with l penalty on the linear coefficients to induce a sparse solution and to set irrelevant parameters to zero after the fitting. this sparse lasso solution is obtained through elastic-net regularized generalized linear models. a supplementary feature selection procedure comparing the full model to a restricted model leaving one tf or mirna out provides the most predominant tf/mirna regulators. tf binding scores are collected from the generic cell line of erythroleukemia cells k from encode for tfs and genes. further, the number of conserved target sites on 3’utr is taken from sequence-based information from targetscan for mirnas and genes. the racer method is implemented in r and publicly available under http://www.cs.utoronto.ca/~yueli/racer.html. the method was evaluated using expression data from an acute myeloid leukemia data set from tcga with samples via a sample based 10-fold cross validation on the prediction of gene expression. to assess the quality of predictions, the spearman rank correlation was calculated resulting in a reassuring value of approximately. further, the full model was compared to models excluding one type of the input variables. the full model performed best and a substantial reduction of spearman correlation was observed by omitting tf regulation and dna methylation. racer also performed with competitive accuracy in predicting known mirna – mrna and tf – gene relationships compared to other methods like genmir++ or encode tf binding scores using e.g., validated interactions from the mirtarbase and knockdown studies. the feature selection procedure revealed predominant transcriptional regulators in the aml dataset. using their associated targets, a functional enrichment analysis showed that dna repair and the tumor necrosis factor pathway were enriched. when applying this panel to cluster patients at different cytogenetic risks, the clustering pattern of the regulatory activities was largely consistent with the risk groups. further, a literature survey on aml showed that many tf regulators among the top predictions had a role in leukemogenesis. rabit regression analysis with background integration is a method for finding expression regulators in cancer by a large scale analysis across diverse cancer types. it integrates tf binding information with tumor profiling data to search for tfs driving tumor-specific gene expression patterns. it can be applied to predict cancer-associated rna-binding protein recognition motifs which are key components in the determination of mirna function.fig. flow chart of rabit method. the input data sets are passed to a linear regression model which calculates sample specific activity values for each regulator and determines general regulatory activities in contrast to our general framework, rabit can, like racer, make use of cnv and dna methylation data additionally integrating promoter cpg content and promoter degree information and takes rbp or tf binding information as regulatory input. the computational model consists of three steps. first, rabit tests in each tumor whether the target genes, identified by the beta method, show differential expression compared to the normal controls including a control for background effects from cnvs, promoter dna methylation, promoter cpg content and promoter degree: gi^=∑fθfbf,i+∑tβtbt,i where gi^ represents the predicted differential gene expression between tumor and normal samples in gene i, b includes values of the f different background factors for gene i, b contains rbp or tf binding information and θ and β are the respective regression parameter vectors. the regression coefficients β are estimated by minimizing the squared difference between measured and predicted gene expression. the regulatory activity score for each tf/rbp is defined by a t-value and its significance by the corresponding t-test. if multiple profiles exist for the same tf from different conditions or cell lines, the profile with the highest absolute value of tf regulatory activity score is selected. in a second step, a stepwise forward selection is applied to find a subset of tfs among those screened in step one optimizing the model error. lastly, tfs with insignificant cross-tumor correlation are removed from the results. computationally, the regression coefficients are calculated via the efficient frisch-waugh-lovell method. tf binding information is taken from tf chip-seq profiles from encode representing tfs and cell types. additionally, recognition motifs for rbps and their putative targets are collected by searching recognition motifs over the 3’utr regions. an implementation of the rabit method can be downloaded from http://rabit.dfci.harvard.edu/download. rabit was applied to tumor profiles of cancer types from tcga using gene expression, somatic mutation, cnv and dna methylation data. to systematically assess the results, the cancer relevance level of a tf was calculated as percentage of tumors with the tf target genes differentially regulated. a comparison to cancer gene databases, i.e., the nci cancer gene index project, the bushman laboratory cancer driver gene list, the cosmic somatic mutation catalog and the ccgd mouse cancer driver genes, showed a consistent picture. further, rabit’s performance was compared to other regression models like lar or lasso where rabit had the best classification results when classifying all tfs into three categories by nci cancer index and achieved better cross-validation error and shorter running time. the regulatory activity of rbps showed that some alternative splicing factors could affect tumor-specific gene expression by binding to target gene 3’utr regions. ismara in contrast to the previous three methods and to our general framework which directly scores tfs or other regulators, ismara infers the activity of regulatory motifs and thereby indirectly deduces the effects of tfs and mirnas. ismara is a web service where no parameter settings or specific processing of the input data, gene expression or chip-seq data are necessary. it can also be used to calculate regulatory activity differences between samples and consider replicates or data from time series.fig. ismara model scheme. the input data sets are passed to a linear regression model which calculates motif activities and determines associated regulators ismara takes sample specific measurements and information about regulatory motifs for tfs and mirnas into account. based on the input of gene expression data or chromatin state measurements, the input signal is calculated for each promoter in each sample. the input signal levels are modelled linearly in terms of the binding site predictions and unknown motif activities: gp,s^=cp+cs+∑mnp,mβm,s where gp,s^ refers to the input signal for a promoter p in sample s, c p and c s are intercepts for each promoter and sample, n p,m summarizes the tf/mirna binding site predictions and β m,s stands for the estimated motif activities. like in the other presented methods, the optimization criterion is to minimize the sum of squared error terms between predicted and measured gene expression. primarily, ismara provides the inferred motif activity profiles sorted by significance and a set of tfs and mirnas that bind to these motifs representing the key regulators. further, a list containing their predicted target promoters, associated transcripts and genes, a network of known interaction between these targets and a list of enriched gene ontology categories is displayed. the web service ismara is available under http://ismara.unibas.ch. ismara employs a bayesian procedure with a gaussian likelihood model and a gaussian prior distribution for β m,sto avoid overfitting. information about regulatory motifs is provided via the annotation of promoters based on deep sequencing data of transcription start sites. to obtain a set of promoters and their associated transcripts, the 5’ ends of mrna mappings from ucsc genome database are clustered with the promoters. tf binding site predictions in the proximal promoter region are collected using position weight matrices representing tfs from jaspar, transfac, motifs from the literature and their own analyses of chip-seq and chip-chip data. additionally, mirna target sites for about seed families are annotated in the 3’utrs of transcripts associated with each promoter. for evaluation, ismara was applied to data from well-studied systems and results were compared to the literature. inferred motif activities were highly reproducible and even more robust than the expression profiles from which motif activities were derived. when comparing samples from human cell types from younger and older donors, ismara was able to identify a key regulator of aging-related changes in expression of lysosomal genes. a joint analysis of the human gnf atlas of tissues and cell lines and the nci- reference cancer cell lines revealed that many of the top dysregulated motifs were well-known in cancer biology like hif1a and has-mir- mirna. they also suggested novel predictions for regulating tfs in innate immunity, mucociliary differentiation and cancer. birte birte takes a mathematically different approach compared to the abovementioned methods integrating tf/mirna target gene predictions with sample specific expression data into a joint probabilistic framework. compared to our general scheme of a tf – gene network, birte takes the tf/mirna – gene network without the interactions between regulators to estimate regulatory activities and infers the network between regulators in a second step.fig. scheme of birte method. the input data sets are passed to a likelihood model which determines active regulators birte takes as input differential gene expression data, an underlying regulatory network including tf/mirna – target gene binding information and optionally cnv data, mirna and tf expression measurements. as opposed to our general framework, birte defines a likelihood model for the set of active tfs/mirnas based on the entire gene expression data d and certain model parameters θ: ld,θr=pd|r,θ=∏d^pd^|r,θ=∏d^∏c∏ipd^ic|rc,θ here d represents the set of all available experimental data including mrna, cnv, mirna and tf expression data and d ic refers to its ith feature measured under experimental condition c. the condition specific hidden state variables r c are estimated with help of the markov chain monte carlo method where a regulator can switch from an active to an inactive state or an inactive and an active regulator exchange their activity states. thereby, the posterior probability for each regulator and condition to influence the expression of its target genes is estimated. simultaneously, a variable selection procedure is applied to achieve sparsity of the model. the optimization goal is not, as one would expect, to return the configuration with highest posterior probability among all sampled ones but to take marginal selection frequencies during sampling into account and filter those above a defined cutoff. after the determination of active regulators, the associated transcriptional network containing tfs and mirnas is inferred from the observable differential expression of target genes and target gene predictions for individual regulators. in practice, the stochastic sampling scheme based on mcmc allows swap operations only when regulators show a significant overlap of regulated targets. the variable selection procedure is implemented via a spike and slab prior which can integrate prior knowledge about the activity of regulators. to infer the associated transcriptional network, nested effects model structure learning is applied. an input mirna – gene network is constructed based on mirmap for mirnas. the tf – target gene network with tfs is compiled by computing tf binding affinities to promoter sequences according to the trap model using data from ensembl, transfac, jaspar and metacore™. an implementation of birte is available for r on bioconductor under https://bioconductor.org/packages/release/bioc/html/birte.html. several simulations were conducted to study model behavior. on the basis of a human regulatory sub-network and accordingly simulated expression data of target genes birte was compared to birta, gemula and a hypergeometric test and further to other network reconstruction algorithms like aracne, genie and genenet. birte performed best in regulator activity predictions including a favorable computation time and was robust against false positive and false negative target gene predictions. additionally, birte was applied to an e.coli growth control and to a prostate cancer data set including normal and cancer samples from geo with corresponding array data from human mirnas and the results showed a principal agreement with the biological literature. aracne we compare aracne as an established, yet local, tool for the reconstruction of gene regulatory networks to the previous five recent genome-wide approaches. the algorithm is background knowledge-free and identifies transcriptional interactions based on mutual information including non-linear and non-monotonic relationships and distinguishes between direct and indirect relationships. aracne is a free tool available under http://califano.c2b columbia.edu/aracne.fig. aracne flow chart. the input data set is used to calculate pairwise mutual information where indirect interactions are removed and which allow a reconstruction of the gene regulatory network aracne uses as input only microarray expression profiles and estimates candidate interactions by calculating the pairwise gene expression profile mutual information i defined as i=ii,j=s+s−s where s denotes the entropy. i i,j measures the relatedness of genes g i and g j and equals zero if both are independent. in a second step, the mutual information values are filtered using an appropriate threshold depending on the distribution of all mutual information values between random permutations of the original data set and indirect interactions are removed. computationally, a gaussian kernel operator is used to calculate mutual information scores. in a subsequent step, the data processing inequality is applied to remove probably indirect candidate interactions. the dpi states that if the genes g i and g k interact only through a third gene g j, then i≤min,i) thus, the least of the three mutual information scores can come from indirect interactions only. aracne’s performance was evaluated on the reconstruction of realistic synthetic datasets and on an expression profile dataset consisting of about b lymphocytes derived from normal, tumor-related and experimentally manipulated populations against relevance networks and bayesian networks. regarding the synthetic networks, aracne had consistently better precision and recall values compared to the two other algorithms and reached very good precision at significant recall levels. it recovers far more true connections and fewer false connections than the other methods with better performance on tree-like topologies compared to scale-free topologies. a reconstructed b-cell specific regulatory network was found to be highly enriched in known c-myc targets where about 50% of the predicted genes to be first neighbors were reported in the literature.  we described five recent methods for the genome-wide inference of regulatory activity, namely the approach by schacht et al., racer, rabit, ismara, and birte. they all assume the topology of the regulatory network to be known, cast activity estimation as an optimization problem regarding the difference between predicted and measured values, take different types of sample specific omics data into account, and eventually produce a list of regulators like transcription factors or mirnas, ranked by their estimated activities in the samples under study. we also included aracne which is background knowledge-free and uses only local dependency measures to reconstruct a regulatory network and indirectly infer activities. all of the presented methods essentially follow the same goal, i.e., accurate ranking of regulatory activity, but differ in the types of measurements being integrated, the background knowledge necessary for their application, the complexity and refinement of the underlying model of gene regulation, and the concrete paradigm used for solving the optimization problem. most of the methods, except for the approach by schacht et al., are available online via a downloadable implementation, a web service, or an r package providing an operable solution for the interested user. whereas an overview of the main features of each method ca be found in table  we now first compare the algorithms regarding their general properties in a descriptive way. the data sets used for evaluation vary between all methods. therefore, we further implemented an evaluation framework to compare the method by schacht et al., racer, rabit and birte in an objective and quantitative way. we used experimental data of three publicly available data sets from tcga and a regulatory network as background knowledge. we first used only mrna expression data as input to the four methods to ensure the result’s comparability, whereas in a second evaluation step, also other omics data sets were included where possible. we further analyzed the relevance of regulators found by different methods using a literature search. general properties experimental data types included the methods differ in the types of measurements being integrated, which corresponds to the level of detail of their model of gene regulations. all six methods use mrna as input. racer, rabit and birte can also integrate cnv, dna methylation, tf/mirna expression data, or somatic mutations. ismara calculates an input signal from microarray, rna-seq, or chip-seq data. additionally, all presented methods use prior knowledge about the underlying regulatory network. these networks are extracted from different data sources and pre-processed in different manners. all methods require at least knowledge about tf – gene relationships, yet racer, birte and ismara also incorporate information about mirnas. when using rabit, the user can choose whether to provide tf or rna binding protein information. the approach of schacht et al. and birte extract regulatory information partly from the commercial metacore™ database, whereas the other methods use only publicly available databases, like encode, jaspar or transfac. the networks which are used for the evaluations published in the respective papers are publicly available for the case of racer, rabit and birte. neither schacht et al. nor ismara make this data available. mathematical models of regulatory activity the methods use different mathematical models to infer regulatory activity. the approach by schacht et al., racer, rabit and ismara use linear regression whereas birte applies a probabilistic framework. aracne, as a local method, is based on mutual information. racer and rabit can be seen as extensions of the approach by schacht et al. since they essentially use the same model structure but incorporate more input data types and more classes of regulatory information. further, racer applies a two-stage regression to infer regulatory activity. optimization frameworks for assessing regulator activities, schacht et al., racer, rabit and ismara minimize the sum of error terms between measured and predicted gene expression. however, the methods use rather different algorithms for solving the resulting optimization problem, and also apply different constraints to achieve model sparsity, robustness of inference, and feature selection. in the approach by schacht et al., the regression model is computed for each gene separately and allows only a maximum number of six regulating tfs. racer uses a lasso approach, while ismara follows a bayesian model that infers regulator activities as posterior distributions. lasso can be interpreted as a bayesian model using laplacian priors instead of gaussian priors in the regression framework obtaining point estimates of the regulatory activities and enforcing sparseness of the solution. in contrast, birte uses a likelihood model with a spike and slab prior to induce model sparsity. this approach implements a selective shrinkage of model coefficients such that estimates are less biased compared to a lasso prior. with the help of the spike and slab prior, sparsity can be controlled in a variable dependent manner allowing the inclusion of prior belief in the activity of each regulator. computed outputs schacht et al. and birte determine activity of regulators over all samples at once, whereas racer and birte first infer sample-specific activities which are combined to cross-tumor activities only in a second optimization step. in contrast, ismara in first place infers motifs activity; these activities are used to deduce the effects of tfs and mirnas by their motif binding profiles. ismara primarily provides sample specific tf and mirna activity but also offers an option to group samples and compare average regulatory activity between different conditions. like birte and aracne, it also infers the network of the regulators themselves. methods and data sets used for evaluation the type and extent of evaluation performed for the different methods vary greatly. they range from direct application to biological problems over the comparison of results to the biological literature to simulation studies. all methods published evaluations results on publicly available datasets, e.g., from the national cancer institute, tcga or geo, but unfortunately address different tissues and cancer types. sample-based cross-validation is applied in the work by schacht et al., racer, rabit and ismara. the first two of these methods use correlation coefficients between measured and predicted gene expression for assessing prediction quality. racer, rabit and birte compare their results to the outcome of other algorithms and to those of restricted models, for example excluding one type of the input variables. all methods search the literature to compare their predictions to previously published studies on the respective biological question. overall, ismara provide the most extensive biological evaluation using a battery of relevant use cases, whereas birte excels in systematic simulation studies. sadly, there are very few works which compare any of the methods presented on the same problem; the only result we are aware of compared aracne and birte regarding their performance in network reconstruction on simulated data, in which birte attained higher robustness against false positive and false negative target gene predictions. quantitative comparison although certain evaluation steps were carried out for all methods, results in the original papers are not comparable as they used different input datasets, different background regulatory networks, and different evaluation metrics. therefore, in addition to the comparison of general properties of the methods, we implemented an evaluation framework using three independent and publicly available test data sets to compare the method by schacht et al., racer, rabit and birte in an objective and quantitative way. all evaluated methods were given the same regulatory network as input. data sets for the evaluation we used experimental data from tcga for three cancer types: colon adenocarcinoma, liver hepatocellular carcinoma and pancreatic adenocarcinoma. for all three cancer types, mrna expression, cnv, dna methylation and mirna expression data is available for primary tumor and normal tissue samples. these data sets are openly accessible via the nci genomic data commons data portal or the nci genomic data commons legacy archive. for mrna gene expression we used processed rna-seq data in the form of fpkm values. the files included ensembl gene ids which were converted to hgnc symbols using the ensembl biomart tool to match the ids of the tf – gene network. in two cases, when multiple ensembl gene ids mapped to one hgnc symbol, we chose the gene with highest log fold change between case and control group. mirna expression was given as rpm measurements. both mrna and mirna data were centered using a weighted mean such that the mean of the case group equaled the negative mean of the control group, and normalized via a weighted standard deviation. cnv data was retrieved as masked copy number segment where the y chromosome and probe sets with frequent germline copy-number variation had already been removed. chromosomal regions were mapped to genes using the r package biomart. if multiple records mapped to one gene, the median of the segment mean values was calculated. for dna methylation data we used the beta-values of illumina human methylation arrays as methylation scores. multiple scores for the same gene were averaged within a sample. we restricted our analyses to the samples for which all four input data types were available. when multiple measurements for one sample and data type were available, we used only the first one in alphabetical order of the file name. after this selection procedure, samples remained for coad, for lihc and for paad. a list including sample and file information is available in additional file together with the experimental data, all evaluated methods were given the same regulatory network as input. we used a publicly available human tf – gene network based on a text-mining approach and complemented it with tf – gene interactions from the public transfac database. this network included interactions between tfs and genes. the network is provided in additional file evaluated methods we conducted the quantitative comparison for the method proposed by schacht et al., racer, rabit and birte. ismara was not included since it is only available as a web service, can only be used with its own, proprietary underlying regulatory network model, and requires the upload of raw data which is prohibited by tcga’s terms of use. also aracne was not included in the quantitative evaluation since it does not use background knowledge and we therefore consider its results as incomparable to the other methods.for the approach by schacht et al. we re-implemented their method as closely as possible to the original design using python and the cuneiform workflow language. due to the high number of integer parameters in the original method, the complexity of optimizing the whole network at once would have by far exceeded computational measures. therefore, like in the original paper, we computed the model for each gene separately and restricted the number of regulating tfs per gene to six. we added a second step where we used these tf – gene interactions building a sub-network to optimize tf activity globally to describe the interplay of the tfs’ effects on their target genes. as in the implementation of schacht et al., we used the gurobi optimizer.7 for racer we used the available r scripts and extracted the resulting sample-specific regulatory activities. rabit published a c++ implementation which they provide on their website and which we used with the fdr option set to as rabit takes differential expression into account, we used the difference of expression values between case and control group as input and ordered the tfs by t-value as proposed in the rabit paper. birte is available as a bioconductor r package. we used r version. with birte version. and applied the method “birtelimma” to estimate regulatory activities with the options niter and nburnin set to as birte has a randomized component, the resulting tf activities are not exactly the same for different runs. we averaged the final activity scores over iterations of birtelimma. for our re-implemented method by schacht et al. and racer we computed separate models for case and control group and ranked the tfs by their activity difference between the two groups. to ensure the result’s comparability, we first used only mrna expression data as input to the four methods. in a second evaluation, we included also other omics data sets where possible. birte was evaluated on mrna and cnv data, rabit on mrna, cnv and dna methylation data, and racer additionally used mirna expression as input. we obtained lists with the regulators ranked according to the absolute value of their computed activity for each cancer type and method, with and without the use of additional inputs. for each cancer type we calculated the size of the overlaps in the four different results using the top and top regulators. the results for the top regulators using either only mrna or multiple omics data sets as input are shown in table  table hgnc symbols of the top regulators found by each method for coad, lihc and paad and the use of only mrna data as input and multiple input data sets. tfs with equal activity values are marked with*. tfs found by several method’s top are marked in bold, blue, red or yellow only mrna as input when only mrna is used as input, one tf is commonly found by the three methods racer, rabit and birte in each data set, respectively: phox2b for coad, epas for lihc and elf for paad. a literature search of these tfs and their targets revealed clear associations to the respective cancer type. the tf obtained commonly for coad, phox2b, is related to tlx a gene which has been shown to play a role in the tumorigenesis of gastrointestinal stromal tumors. epas which was found in the lihc top tfs of three methods, is linked to cxcl which plays an important role in metastasis formation of hepatocellular carcinoma by promoting the migration of tumor cells. for paad, three methods ranked tf elf high, which is related to genes in our network, inter alia to brca and lyn. mutations in the brca gene have been implicated in pancreatic cancer susceptibility, whereas the knockdown of lyn reduced human pancreatic cancer cell proliferation, migration, and invasion. these results underline that the methods are able to find biologically relevant information about regulation processes in cancer. several tfs in the top are found by two of the four methods for instance, racer and rabit have four common top tfs in the coad data set. however, the top tfs found by the method by schacht et al. do not overlap with any top tfs of the other methods in any data set. the agreement of racer, rabit and birte in the top tfs hints to the biological importance of the found tfs since this overlap is statistical significant as the probability of finding common tfs in three sets of ten randomly chosen ones out of tfs is below. additionally, the methods do identify different tfs for different data sets, indicating the importance of the actual cancer specific mrna expression values and that results are not dictated by the background network. the results for the number of overlapping regulators in the top between the four methods and the three different data sets are shown in fig.  for rabit, only tfs for coad could be ranked since all other tfs had an activity value equal to zero.fig. number of overlapping tfs in the top of ranked tfs per method in coad/lihc/paad is shown) when looking at the overlap of three of the four methods, the number of overlapping tfs is still the highest for the triplet racer, rabit and birte. for the lihc dataset two tfs are found in the top of all four methods. e2f is a downstream target of zbtb which was associated to the expression of cell cycle-associated genes in liver cancer cells. two target genes of e2f cdk and tp were also involved in liver cancer development and proposed as prognostic marker of poor patient survival prognosis in hepatocellular carcinoma. further, epigenetic alterations of the ednrb gene, a target of sox might play an important role in the pathogenesis of hepatocellular carcinoma. even if the result of four methods finding two common tfs is not statistically significant, their association to liver hepatocellular carcinoma shows that the methods reach their goal of identifying relevant tfs. however, when comparing different data sets, the methods tend to rank the same tfs under the top to a greater or lesser extent. for example, the overlap of all top tfs of the three cancer types is only one tf for rabit and nine tfs for birte, but tfs for the method by schacht et al. and even tfs for racer. therefore, the results from rabit and birte seem to be more cancer type specific and less dependent on the regulatory network than the results from racer. however, we did not specifically investigate the influence of the underlying network and its topology on the results which would be an interesting point for further research. multi-omics data as input when not only taking mrna into account but also mirna, cnv and dna methylation, the results are more difficult to compare between the methods, since they all use a different way of combining different types of data due to their models and implementations. we are aware of the lower level of comparability of this approach regarding the multi-omics results in contrast to a scenario, where all methods are evaluated on the same set of input data. however, we intended to use maximum set of input data for each method to cover the effect of the use of multiple omics data sets compared to only mrna as input. birte was evaluated on mrna and cnv data, rabit on mrna, cnv and dna methylation data, and racer additionally used mirna expression as input. whereas racer and rabit considered cnv or dna methylation data as one background factor and compute only one activity value, birte evaluated the influence of each cnv separately. the results show that racer exclusively ranks mirnas high; not a single tf is found among the top regulators. also, the influence of cnvs was high in lihc and paad. however, the tfs that racer found in the top when using only mrna data as input are still ranked high in the multi-omics scenario, e. g the coad top three tfs of the mrna results are ranked 13th, 16th and 14th in the results of the multi-omics input. the difference of the results coming from the two input types is less for rabit: seven tfs are still in the top for coad when using cnv and dna methylation additionally to mrna data. therefore, the contribution of additional input data seems not to be crucial for the performance of rabit. birte considers each cnv as a potential regulator which increases the total number of regulators enormously. still, two commonly present tfs in the top of the coad data set are found by either the sole mrna input and the multi-omics approach. the overlap of the top of rabit and birte in the multi omics case is considerable with three tfs in lihc, and one tf in paad. three of them were already found when using only mrna data as input. the results for the use of different input data sets show that the top ranked regulators are drastically changed when using additionally mirna data in racer, but change less when only cnv or dna methylation data is provided in rabit and birte. however, the results from multi omics analyses are difficult to compare since the combination of input data sets is not consistent across the three different methods. discussion background networks a crucial input to the models is the underlying regulatory network which is needed to reduce the search space for actual regulatory activity. however, the construction of comprehensive tf/mirna – gene regulatory networks is difficult for various reasons. firstly, a comprehensive characterization of the human regulatory repertoire is lacking since only about half of the estimated – tfs in the mammalian genome is known. chip experiments, prone to a high false positive rate, were used to identify tf binding patterns but each assay is limited to the detection of one tf in one condition and therefore tf binding has not been characterized for many tfs in most cell types. further, the local proximity of a binding site to the transcriptional start site of a gene does not automatically implicate transcriptional regulation. with regard to posttranscriptional regulators, the functions for only a few of the around different mirnas have been experimentally determined and current data on mirna targets is mostly based on computational predictions. generally, the knowledge about tf and mirna binding is scattered over the biological literature and different, partly commercial, databases, impeding the construction of comprehensive networks. therefore, any comparative evaluation of the methods presented here would have to make sure that the same background network is used for each computation. besides, studies on the impact of network incompleteness or different error rates in networks would be important to assess the ability of the methods to cope with such common problems. simulation studies will be vital in this regard. the graph view on regulation the modelling of regulatory networks as graphs, as used in all presented methods, is perhaps not the optimal representation for the underlying biological regulatory processes. a graph cannot easily account for important effects such as tf complex formation and temporal and spatial synchronization of activities. furthermore, tf binding is affected by chromatin state and the impact of posttranslational modifications on transcriptional activity which are difficult to include in a graph view on regulation. the model’s dependence on the topological structure and the robustness to changes in the underlying network have not been evaluated or discussed in any of the presented methods even if these issues are known to have a severe influence in network analysis. underlying mathematical model linear models, widely spread in different fields of science, provide a simple and easily understandable design but over-simplify the underlying biological processes. nonlinear behavior, e.g., saturation effects, cannot be represented. considering that the number of available samples is typically relatively small, the incorporation of many different data types and according parameters into the model could result in excessively complex designs prone to overfitting, but this issue lacks general awareness. only two of the presented methods incorporate parameter priors, and two apply cross validation techniques to estimate prediction performance. further, the effect of temporal buffering between tf binding and the actual effect on gene expression is not included in any of the methods. comparability all methods produce a ranked list of regulators. comparing these results across different methods, even when applied on the same data set and using the same background network, is difficult since no generally accepted benchmarks are available. therefore, there currently is no objective measure to designate a best method. the closest comparable evaluation effort we are aware of is implemented in the “dream – network inference” challenge, which targets gene regulatory network reconstruction. the invited participants reverse-engineered a network from gene expression data, including a simulated network, and evaluated the results on a subset of known interactions or the known network for the in-silico case. the approach of genie which trains a random forest to predict target gene expression performed best and the integration of predictions from multiple inference methods showed robust and high performance across diverse data sets. however, an extensive competitive evaluation to determine active regulators based on a given regulatory network has, to the best of our knowledge not been carried out yet. we therefore compared the results of four methods in a quantitative way. the experimental data and the regulatory network we used as input are publicly available to ensure transparency of our results. the results suggest that the methods are able to find biologically relevant information about regulation processes in cancer. however, the result overlaps are rather low. this seems surprising as all methods essentially follow the same goal, i.e., identification of the most differentially active tfs or genes. we think further research is necessary to exactly characterize the specific strengths of each method. furthermore, we did not investigate the influence of the underlying network on the results, which is another topic for further research.  despite their often rather involved procedures and models, none of the presented methods adequately reflects the biological reality of regulatory activity in cells. a specific disease phenotype is rarely caused by a single gene but rather a product of the interplay of genetic variability, epigenetic modifications and post-transcriptional regulation mechanisms. the presented methods ignore a multitude of such factors like the effects of chromatin state and alternative splicing, nonlinear relationships between regulatory activity and gene expression, or kinetic and temporal effects. furthermore, tfs themselves regulate the expression of other tfs forming feedback loops which are not considered in any of the presented methods. nevertheless, the methods apparently are able to detect strong signals and produced promising results in terms of ranking transcription factors by their activity and are thus valuable tools for identifying biomarkers for specific phenotypes. additional files additional file 1: lists information about the samples and files from tcga included in our quantitative evaluation for all three cancer types. additional file 2: includes an adjacency list of the connected nodes of the tf – gene network. the list includes three columns where each row indicates an association with the value of “edge” between a tf and a gene. complexes of tfs are indicated with a separating “.” between their components. http://www.gurobi.com/products/gurobi-optimizer http://lsresearch.thomsonreuters.com/pages/solutions/1/metacore https://gdc-portal.nci.nih.gov https://gdc-portal.nci.nih.gov/legacy-archive http://www.ensembl.org/biomart/martview, release 87 http://www.gene-regulation.com/pub/databases.html, release version, available under a free academic license http://www.cs.utoronto.ca/~yueli/racer.html http://rabit.dfci.harvard.edu 
furthermore, we quantitatively compare the results of four of the presented methods based on publicly available data sets. however, despite their common root in mathematical optimization, they vastly differ in the types of experimental data being integrated, the background knowledge necessary for their application, the granularity of their regulatory model, the concrete paradigm used for solving the optimization problem and the data sets used for evaluation. in various settings, these methods produced promising results in terms of estimating transcription factor activity and identifying key biomarkers for specific phenotypes. the results show that all methods seem to find biologically relevant information. here, we review five recent methods of this class in detail and compare them with respect to several key properties. recently, several algorithms have been proposed which model genome-wide gene regulation as sets of equations over the activity and relationships of transcription factors, genes and other factors. gene regulationregulatory networksystems biologymathematical optimizationhttp://dx.doi.org/ /501100001659deutsche forschungsgemeinschaft10911091trescher saskia münchmeyer jannes issue-copyright-statement© the author 2017 gene regulation is one of the most important cellular processes, indispensable for the adaptability of organisms and closely interlinked with several classes of pathogenesis and their progression.
 gene regulation is one of the most important cellular processes, indispensable for the adaptability of organisms and closely interlinked with several classes of pathogenesis and their progression. elucidation of regulatory mechanisms can be approached by a multitude of experimental methods, yet integration of the resulting heterogeneous, large, and noisy data sets into comprehensive and tissue or disease-specific cellular models requires rigorous computational methods. recently, several algorithms have been proposed which model genome-wide gene regulation as sets of equations over the activity and relationships of transcription factors, genes and other factors. subsequent optimization finds those parameters that minimize the divergence of predicted and measured expression intensities. in various settings, these methods produced promising results in terms of estimating transcription factor activity and identifying key biomarkers for specific phenotypes. however, despite their common root in mathematical optimization, they vastly differ in the types of experimental data being integrated, the background knowledge necessary for their application, the granularity of their regulatory model, the concrete paradigm used for solving the optimization problem and the data sets used for evaluation.  here, we review five recent methods of this class in detail and compare them with respect to several key properties. furthermore, we quantitatively compare the results of four of the presented methods based on publicly available data sets.  the results show that all methods seem to find biologically relevant information. however, we also observe that the mutual result overlaps are very low, which contradicts biological intuition. our aim is to raise further awareness of the power of these methods, yet also to identify common shortcomings and necessary extensions enabling focused research on the critical points. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords gene regulationregulatory networksystems biologymathematical optimizationhttp://dx.doi.org/ /501100001659deutsche forschungsgemeinschaft10911091trescher saskia münchmeyer jannes issue-copyright-statement© the author 2017 