transfactor: Transcription factor activity estimation via probabilistic gene expression deconvolution

transfactor: Transcription factor activity estimation via probabilistic gene expression deconvolution | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results transfactor: Transcription factor activity estimation via probabilistic gene expression deconvolution View ORCID Profile Koen Van den Berge , Peter Bickel , View ORCID Profile Sandrine Dudoit doi: https://doi.org/10.1101/2025.03.19.644088 Koen Van den Berge 1 Department of Statistics, University of California , Berkeley, CA, USA 2 Department of Applied Mathematics, Computer Science and Statistics, Ghent University , Belgium Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Koen Van den Berge Peter Bickel 1 Department of Statistics, University of California , Berkeley, CA, USA 3 Graduate Group in Biostatistics, University of California , Berkeley, CA, USA 5 Center for Computational Biology, University of California , Berkeley, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sandrine Dudoit 1 Department of Statistics, University of California , Berkeley, CA, USA 4 Division of Biostatistics, School of Public Health, University of California , Berkeley, CA, USA 5 Center for Computational Biology, University of California , Berkeley, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sandrine Dudoit For correspondence: sandrine{at}stat.berkeley.edu Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Gene expression is the primary modality being studied to differentiate between biological cells. Contemporary single-cell studies simultaneously measure genome-wide transcription levels for thousands of individual cells in a single experiment. While the characterization of cell population differences has often occurred through differential expression analysis, tiny effect sizes become statistically significant when thousands of cells are available for each population, compromising biological interpretation. Moreover, these large studies have spurred the development of methods to infer gene regulatory networks (GRNs) directly from the data, and GRN databases are becoming more comprehensive. In this work, we propose a statistical model for gene expression measures and an inference method that leverage GRNs to deconvolve transcription factor (TF) activity from gene expression, by probabilistically assigning mRNA molecules to TFs. This shifts the paradigm from investigating gene expression differences to regulatory differences at the level of TF activity, aiding interpretation and allowing prioritization of a limited number of TFs responsible for significant contributions to the observed gene expression differences. The inferred TF activities result in intuitive prioritization of TFs in terms of the (difference in) estimated number of molecules they produce, in contrast to other widely-used methods relying on arbitrary enrichment scores. Our model allows the incorporation of prior information on the regulatory potential between each TF and target gene through prior distributions, and is able to deal with both repressing and activating interactions. We compare our approach to other TF activity estimation methods using two simulation experiments and two case studies. 1 Introduction In single-cell RNA sequencing (scRNA-seq), researchers investigate gene expression levels as driving factors of observed biological differences between cell types, disease states, and responses to treatments. The resolution of gene expression at the level of single cells has been a boon for biological research, as it has allowed to deconvolve cellular gene expression heterogeneity, as well as answer a new set of questions, including the identification of homogeneous cell entities (‘cell types’) or dynamic cellular differentiation processes (‘trajectories’) such as cell development and activation. Traditionally, differences between cell types or within/between trajectories have been identified using differential expression (DE) analysis, where gene-level differences in average expression are assessed using hypothesis testing [ Finak et al., 2015 , Van den Berge et al., 2020 ]. Through this approach, the study and interpretation of differentially expressed genes (DEG) therefore provides an opening to elucidate the biological processes of interest. However, the large number of sequenced cells (in the tens or hundreds of thousands) of contemporary scRNA-seq datasets renders biologically irrelevant (or insignificant) effect sizes statistically significant. While moderation methods such as those implemented in the limma framework [ Smyth, 2004 ] were designed in part to help with the issue of biologically uninteresting statistical significance, they are unable to fully resolve the problem. This hampers interpretation and leads to ad hoc tuning of significance thresholds as low as 10 — 12 [ McFaline-Figueroa et al., 2019 ]. Even then, the fraction of genes that are in the set of DEG is often too large for reasonable interpretation by biologists. One possible solution is to test against a fold-change threshold [ McCarthy and Smyth, 2009a ], but this again relies on the ad hoc selection of an appropriate threshold for each hypothesis being assessed. Gene expression is a molecular process that is being regulated by transcription factors (TFs), that is, proteins that bind DNA regions to be transcribed, resulting in the production of pre-messenger RNA (pre-mRNA) molecules, in the case of an activating regulator, or the inhibition of the production (e.g., through spheric hindrance) of premRNA molecules, in the case of a repressing regulator. A more parsimonious way to interpret gene expression differences is to move one level up in the hierarchy: instead of identifying which genes are differentially expressed, identify transcription factors that are driving the observed gene expression differences. Since there are far fewer TFs than genes, this leads to a more manageable interpretation in the hands of biologists. For example, a recent in-depth review by [ Lambert et al.2018 ] provides a manually-curated list of human TFs and identified 1,639 proteins to be “known or likely TFs”, while the current human genome annotation corresponds to around 20,000 protein-coding genes [ Frankish et al., 2019 ]. However, while TF protein abundance is typically high in single cells, the RNA concentration of the respective TF genes is often low [ Weinreb et al., 2020 , Chen et al., 2021 ], and the activity of TF proteins may be modulated through biochemical processes such as phosphorylation [?]. It is thus possible that TFs that are highly active, i.e., producing many RNA molecules from their downstream target genes, have relatively low RNA concentration. This effectively jeopardizes direct TF activity estimation using the corresponding TF’s RNA transcripts, motivating statistical methods to estimate TF activity using downstream target genes instead. Several methods have previously used (sc)RNA-seq data for identifying TFs responsible for driving gene expression differences, while others focus on alternative data modalities such as single-cell chromatin conformation data [ Berest et al., 2019 ]. Gene-expression-based methods typically either first estimate a gene regulatory network (GRN) or require one as input. A GRN specifies the genes being (positively or negatively) regulated by each TF. A set of genes being regulated by one TF is also often called a regulon . Some methods, such as TFEA [ Rubin et al., 2020 ], start from a list of regions of interest based on, for example, chromatin conformation or ChIP-seq data, and perform an enrichment analysis to assess which motifs are more abundant in the regions of interest than what could be expected by chance. Other methods, such as viper [ Alvarez et al., 2016 ] and AUCell [ Aibar et al., 2017 ], use a gene regulatory network to estimate TF activity for each TF directly using the gene expression data. For all three methods, however, the returned activity scores are on arbitrary scales and are therefore not directly interpretable. In this manuscript, we introduce a new method that uses scRNA-seq data to infer TF activity in terms of readily interpretable estimates of the number of mRNA molecules produced by each TF for each gene in a particular cell. The method, called transfactor, requires a matrix of gene expression measures and a gene regulatory network as input, and probabilistically deconvolves overall gene expression measures into TF-specific gene expression measures that reflect the allocation of transcripts to TFs that may have produced them. transfactor is capable of handling both activating as well as repressing edges in a GRN. A Bayesian extension of the model allows to incorporate prior knowledge on the confidence of each edge in the network, while also allowing the estimation of these prior parameters through an empirical Bayes framework. If many edges are expected to be irrelevant for the dataset being studied (for example, when a GRN from a global database is used, not all edges may be relevant for the current experiment), a sparse initialization step of the fitting algorithm ensures sparse results in terms of TF activity measures, effectively boosting TF activity estimation accuracy in these scenarios. We compare several variants of transfactor to state-of-the-art methods, as identified in Holland et al. [2020] , using two simulation experiments, and provide two extensive case study analyses showing how TF activity estimation can unlock biological insight. 2. Methods 2.1 Statistical model 2.1.1 Notation and definitions We assume that we know a gene regulatory network that is relevant for the system under study, which comprises G genes and T transcription factors. We define the G × T adjacency matrix A of the GRN, with elements A gt = 1 ( A gt = —1) if gene g is positively (negatively) regulated by TF t and zero otherwise. In Section 2.2.2 , we show how to account for repressions (negative A gt ) prior to model fitting; in the remainder of the manuscript, we will assume there are no repressing edges in A, i.e., only non-negative entries. Based on A, we know that (i) a gene g is regulated by a known set of TFs which is denoted by 𝒯 g = { t : A gt = 1}, and (ii) a TF t regulates a known set of genes that is denoted by 𝒢 t = { g : A gt = 1}. We assume that cells have been grouped into relevant entities, be it cell classes, states, or experimental conditions, generally corresponding to their cell identity, which we refer to using the general moniker of ‘cell type’. Let n denote the total number of observed cells and n c the number of cells observed for cell type c, c = 1,…, C . For a particular cell i that belongs to cell type c ( i ), let Z gti denote a measure of the number of transcripts of gene g produced by TF t and let denote a measure of the total number of transcripts of gene g . While Y gi is observed in scRNA-seq, Z gti is not, and the goal of our method is to deconvolve Z gti from Y gi . Let µ gtc and µ g . c denote, respectively, the expected number of transcripts of gene g produced by TF t in cell type c and the expected total number of transcripts of gene g in cell type c . Note that The parameter of interest is π gtc , the proportion of transcripts of gene g produced by TF t in cell type c , i.e., where for each gene g and cell type c . We further introduce as the global (i.e., across all g ∈ 𝒢 t ) activity of a TF t in cell type c , defined as the average number of molecules the TF produces across all the genes it is regulating, i.e., Distance-based ranking of transcription factors In order to identify which TFs may contribute most to differences in expression between cell types, we can construct a distance-based ranking of transcription factors based on their activity. Here, we focus on the Euclidean distance, but note that other distance measures could similarly be used. A measure of the contribution of a particular TF t to differences in gene expression between cell types c and d can be defined as To adjust for gene expression magnitude, we can also work with a scaled distance measure Either measure can then be used to rank TFs according to their contribution to differences in expression between cell types. If a comparison across multiple groups is of interest (e.g., in the olfactory epithelium case study), one could sum D t across all relevant comparisons or look at contrasts of the µ gtc . Note that D t is a more informative measure of the contribution of a particular TF t as compared to the Euclidean distance on directly. Indeed, D t also accounts for cases where , but where π gtc ≠ π gtd , i.e., the TF is similarly active overall across cell types, but it is preferentially regulating different sets of genes in the different cell types. 2.1.2 Basic model assumptions For the model as introduced in Section 2.1.1 , the parameters µ gtc are unidentifiable based on only the observed data Y gi . Here, we make two assumptions that should aid towards identifiability. First, we assume a multiplicative model for µ gtc , parameterized according to a gene-specific effect α gc and a TF-specific effect β tc , i.e., Biologically, (A1) means that, for a given gene g , we make the strong assumption that all TFs t ∈ 𝒯 g have the same binding capacity and therefore the competition between different TFs to produce mRNA molecules is only affected by β tc . This is a simplification of reality, where other determinants, such as chromatin conformation, histone marks, and the presence of other TFs, would also affect the TF-specific contribution at a genomic locus [ Bannister and Kouzarides, 2011 , Sönmezer et al., 2021 ]. Further, we assume that This assumption also helps in making the model parameters interpretable, since we can consider the TF-specific contribution in (A1) to correspond to the global activity of that TF. 2.1.3 Implications of model assumptions Our model assumptions, (A1) and (A2), have the following implications that are useful to note explicitly, since they aid in interpreting the model as well as understanding its constraints. First, we note that Equations (3) and (A1) imply that and thus, given (A2), the following constraint must hold From Equations (1) , (A1) , and (A2) , it follows that which allows us to rewrite (A1) in a more readily interpretable and intuitive form which shows that the contribution of a TF t to the gene expression of gene g is related to its activity in that cell type as compared to the activity of all other TFs regulating the gene . This also allows us to relate π gtc to the TF activity , as Equations (2) and (8) imply that i.e., the proportion of transcripts of gene g produced by TF t is only dependent on the TF’s global activity , relative to the activity of all other TFs regulating the gene. 2.1.4 Model identifiability We have not established model identifiability and the proof is a non-trivial mathematical question that is beyond the scope of this paper. The identifiability problem is shared by other methods that analogously attempt to assign k-mers of sequencing reads to transcripts [2016, 2017] (instead of mRNA molecules to TFs). In practice, using small network examples (Supplementary Methods Section 6), we found a few conditions that are necessary, but not sufficient, for model identifiability: T ≤ G . The number of transcription factors must not exceed the number of genes. This condition should always be met, as it is a biological reality. A must be of full column rank, i.e., column rank T . This condition may or may not hold, depending on the provided GRN. If a limited set of TFs are linearly related in terms of their regulated genes (i.e., the corresponding columns in the adjacency matrix are not linearly independent), one may consider removing these TFs in order to obtain a full-column-rank matrix. 2.1.5 Poisson model and a Bayesian extension of the model The basic assumptions (A1) and (A2) of Section 2.1 naturally lead to a method-of-moments estimation approach, with estimating equations which have no closed-form solution for the the parameters α gc and β tc . We will not directly use method-of-moments estimation and will instead introduce additional parametric assumptions that allow us to develop a Bayesian extension of the model that incorporates prior biological knowledge and induces sparsity in the parameter estimates. We assume that the latent variables follow a Poisson distribution, i.e., where the Z gt i are independently distributed across genes, TFs, and cells. This frequentist model will be referred to as the ‘Poisson model’. Note, however, that it does not allow incorporating a priori known information on the edges present in the GRN. Indeed, different edges might have different confidence scores based on empirical evidence supporting each edge, e.g., as obtained via databases or experiments. To incorporate such information, we specify, for each gene g and cell type c , a Dirichlet prior distribution for ( π gtc : t = 1,…, T ), subject to the constraint of Equation (9) , i.e., The ( π gtc : t = 1,…, T ) are assumed to be independently distributed over genes g and cell types c and the prior distributions are constant across cell types. Further, if A gt = 0, then γ gt = 0. This formulation also implies that if genes g 1 and g 2 are regulated by the same TF t , then , and if they are regulated by exactly the same set of TFs ∀ t ), then . The constraint on the Dirichlet distribution that the π gtc have this specific form makes it explicit that there is no inherent sampling variability associated with the Dirichlet prior, but that it is mainly there to borrow information across cell types and/or to incorporate prior knowledge. In our model, the γ gt can be assumed to be either: (i) equal to zero, corresponding to a non-informative (improper) prior distribution, which reduces to inference based on the frequentist Poisson model; (ii) known (e.g., by using information from biological databases on the regulatory strength between a gene and a TF), in which case we are working with a Bayesian model, or (iii) estimated from the data in an empirical Bayes fashion. The model is thus where it follows that Y gi ∼ Poi ( µ g . c ( i ) ) (due to independence over t ) and also . The parameters of interest are π gtc . Note that while downstream interpretation may be focused on either π gtc or , both are readily available once we have estimated one or the other, see Equation (9) . We will also estimate γ gt , if they are not assumed to be known. 2.2 EM algorithm for fitting the hierarchical Poisson model We first note that the parameter µ g . c can be inferred solely from the observed data Y = { Y gi : i = 1,…, n ; g = 1,…, G} and trivially estimated using the sufficient statistic , with C c = { i : c ( i ) = c} , by We then use the expectation-maximization (EM) algorithm [ Dempster et al., 1977 ] to estimate the parameters θ = { π gtc , γ gt }. In the frequentist setting, the algorithm is used to maximize the likelihood function with respect to θ (which, in this case, does not include γ ), and in the (empirical) Bayes setting, it is used to maximize the posterior density with respect to θ (which possibly includes γ as unknown). The frequentist EM algorithm maximizes the Q -function with respect to θ new , where Z = { Z gti : i = 1,…, n ; g = 1,…, G ; t = 1,…, T} . Note that the second term only involves the parameters µ g . c , that have already been inferred from the observed data. Hence, the EM algorithm is only concerned with the first term. In a Bayesian setting, the EM algorithm can be adapted to estimate the mode of the posterior distribution of θ by updating the Q-function in the M-step. Specifically, the EM algorithm proceeds as follows. E-step: Calculate Q by deriving the expected value of Z gti , given the data Y and current estimates θ old . M-step: In a Bayesian setting, maximize Q ( θ new | θ old )+ G ( θ new ) with respect to θ new , where G is the log-prior density [ Dempster et al., 1977 ]. In our case, G corresponds to the log-prior density of the π gtc parameters, which is assumed to be a Dirichlet density with parameters γ gt . Known prior . Maximize the function Q ( θ new | θ old ) + G ( θ new ) with respect to π gtc , keeping γ gt fixed. This corresponds to calculating . Unknown prior . In an empirical Bayes setting, where we also estimate the parameters of the prior distribution, we additionally maximize the function Q ( θ new | θ old )+ G ( θ new ) with respect to γ gt , keeping π gtc fixed at its new value. There is no closed form solution for this, but we will use an approximate MLE as update step. This corresponds to calculating . After each iteration, we set θ old = θ new . Full details on the EM algorithm and update steps are provided in Supplementary Methods. Note that, if the prior parameters γ gt are assumed known, then setting all γ gt = 1 is equivalent to providing no prior information, see Equation (20) in Supplementary Methods. In general, it is not guaranteed that the EM algorithm will converge to the global maximum of the likelihood function, and it is possible that the algorithm will instead converge to a local maximum close to its initialization value. However, initializing the EM algorithm with different random starting values in our simulation studies leads to the same parameter estimates upon convergence, providing empirical evidence of good behavior. We recommend that users investigate diagnostic plots to check convergence of the EM algorithm and, in particular, track the log-likelihood evolution across iterations. Our software package transfactor provides functionality for assessing convergence of the EM algorithm ( Section 2.3 ). 2.2.1 Initialization of the EM algorithm The default initialization for π gtc is , i.e., an equal partitioning of a gene’s transcripts among its regulating TFs. However, while a GRN may indicate that a gene is possibly regulated by many TFs, it is unlikely that all of these plausible TFs are actively contributing to gene expression in a specific biological tissue or condition. It would therefore make sense to obtain sparse estimates of TF activity, where some TFs do not contribute at all to the expression of some genes. If there is prior information, one may promote such sparsity through the Dirichlet prior; indeed, letting γ gt = 0 will enforce TF t to not contribute to the gene expression of gene g . However, this is enforced across all cell types, which could be suboptimal. In addition, prior information to promote such sparsity may not be available. An alternative is to promote sparsity in the TF activity estimates in an unsupervised manner. There are several approaches to do so. A natural paradigm would be to add a sparsity-inducing penalty to the log-likelihood, to be incorporated in the fitting using the EM algorithm. However, many of these approaches result in non-closed-form solutions in the M-step, often requiring numerical optimization at each iteration of the EM algorithm, thereby drastically increasing the computational burden of model fitting. Instead, we use a sparse initialization of the parameters in the EM algorithm. Indeed, if a particular is initialized as , then the update steps of the EM algorithm will keep it zero provided γ gt = 1 (i.e., if we do not feed prior information to the EM algorithm through g, see Equation (20) in Supplementary Methods) or if we use the Poisson model without the Dirichlet prior. In order to derive these sparse initial estimates, we rely on a Gaussian approximation to the Poisson model, as detailed below. A Gaussian approximation to the Poisson model to promote unsupervised sparse solutions The sufficient statistics for µ g. c are Y g. c . We rely on a Gaussian approximation to the Poisson model, where Note that, given Equation (9) and , where, as before, we are interested in estimating (and, correspondingly, π gtc ). The minimum modified Chi-squared objective function to estimate , according to Neyman [1949] , corresponds to where we use a heuristic pseudocount of 1 in the denominator to avoid division by zero. In order to provide sparse initial estimates, we further use a LASSO approach [ Tibshirani, 1996 ], minimizing the penalized empirical risk for each cell type separately, subject to . In the loss function, λ c is a penalty parameter that controls the amount of shrinkage and that is estimated using cross-validation for each cell type separately. The R package glmnet is used for fitting the regression model using the LASSO [ Friedman et al., 2010 ]. The solutions for each cell type will be our initial estimates which will be provided to the EM algorithm. Remark 1 Note that the solutions will be sparse, however, differently so for each cell type, meaning that some TFs will contribute to gene expression in some cell types but not in others. Remark 2 If no prior information is provided, the LASSO initialization above will be decisive in terms of which TFs do not contribute to gene expression within a particular cell type. Indeed, if , then . Remark 3 If prior information is provided, i.e., if γ gt > 1, then the sparsity of the LASSO will be ‘overridden’ by the prior information, see Equation (20) in Supplementary Methods. 2.2.2 Accounting for repressions Our method attempts to assign mRNA molecules to the TFs that produced them. If the GRN contains repressing edges between a TF and a target gene that itself is not a TF, then we discard these interactions (i.e., set A gt = 0). These repressing edges may have hampered the production of molecules for the specific target gene, but this is not something that we attempt to model. Instead, we are only interested in deconvolving counts for molecules that have actually been produced. In contrast, we do account for repressing interactions between TFs, as these may be informative to assign produced molecules to them. Indeed, if a TF’s repressor is highly active, it is less likely that this particular TF will have produced molecules (i.e., it will have a lower TF activity ). In order to account for this type of repression, we apply the following procedure prior to fitting the (hierarchical) Poisson model. A TF that is itself being repressed by other TFs will occur both as a row and as a column in the adjacency matrix A . We denote the row of A corresponding to the gene g t for TF t by and the set of TFs that are repressing TF t by . Also, let denote the (observed) gene expression measure of repressed TF t in cell i from cell type c ( i ) and denote the sum of the gene expression measures across all TFs that are repressing TF t . If, in the system under study, TF t is actually repressed by (a subset of) the TFs in R t , we expect to be negatively associated with . We model as a function of log using a Poisson generalized linear model In this model, a negative β 1 corresponds to repressing interactions. We choose to account for repressing interactions if β 1 is found to be significantly different from zero at nominal significance level of 5% after fitting the model in Equation (17) . Repressing interactions are then handled by multiplying using a scaling factor ρ tc ∈ [0, 1] estimated as reflecting the ratio of the expression of TF t given the average expression of its repressing TFs relative to no expression of its repressing TFs. Here, . See Supplementary Methods for details on how we use ρ tc in the model fitting. 2.3 Software implementation The proposed methodology is implemented in the open-source R software package transfactor, which will be submitted to the Bioconductor Project. TF activity can be estimated using the (hierarchical) Poisson model with or without information for the Dirichlet prior distribution. In the fitting of the EM algorithm, we work with sufficient statistics Y g. c to ease the burden on computational memory requirements. 2.4 Benchmarking and case studies 2.4.1 Evaluated methods We evaluate three different variants of our method: poisson: Poisson model. dirMult: Hierarchical Poisson model with known Dirichlet prior. dirMultEBayes: Hierarchical Poisson model with unknown Dirichlet prior. For each of these variants, the EM algorithm is initialized either with or without the sparse LASSO-based approach described in Section 2.2.1 , which we indicate by the suffix lasso , leading to a total of six different methods to evaluate. Except for the Gamma-model-based simulation framework from Section 2.4.3 , where no repressing interactions are present, we also evaluate the effect of accounting for repressions as described in Section 2.2.2 , which we indicate by the suffix repr . As competing methods, we also evaluate viper [ Alvarez et al., 2016 ] and AUCell [ Aibar et al., 2017 ], which are, top-performing methods according to the benchmarking study of and Holland et al. [2020] . viper . We use viper v1.24.0. We run a single-sample analysis, where TF activity is estimated for each sample (cell) independently, by calculating normalized enrichment scores using the scale method. In brief, viper compares the gene expression of each sample to the average across all samples. The expression counts for each gene are standardized to have zero mean and unit standard deviation across cells. A rank-based enrichment score is computed, where it is assumed that a TF is active if its targets have high ranks (i.e., are highly expressed). The algorithm is capable of incorporating both activating as well as repressing edges, and also information on the confidence of each edge between a TF and a target gene. We run viper using the same confidence scores as the dirMult model of transfactor, effectively giving it an advantage as compared to other methods not using these confidence scores. AUCell . We use AUCell v1.12.0. Like viper , AUCell also first scales the expression counts to a gene-wise mean of zero and standard deviation of one across all cells/samples before calculating TF activity scores. AUCell calculates TF activity scores by ranking the (by default, top 5%) scaled gene expression measures within each cell, and calculating the area under the receiver operating characteristic (ROC) curve for each regulon (i.e., set of genes regulated by a TF), where a ‘true positive’ is considered a gene that belongs to the regulon. The enrichment score thus reflects the presence of genes from the regulon in the top-expressed genes in the cell. The main difference between transfactor and competing methods AUCell and viper , is that transfactor is based on a statistical model that aligns with the data-generating mechanism and makes clear what the assumptions are. Through this, the parameter estimates of transfactor are interpretable; for each gene, one can investigate the relative contribution of TFs to its expression, and how this possibly varies across conditions. 2.4.2 Differential activity analysis In order to identify TFs that are differentially active under different conditions, we first obtain the estimated number of molecules produced by each TF (Equation (21) in Supplementary Methods) using transfactor . We then model the estimated counts using a negative binomial generalized linear model as implemented in glmGamPoi [ Ahlmann-Eltze and Huber, 2020 ]. The design matrix corresponds to the conditions of interest, i.e., cell type groups or pseudotime bins, and inference is carried out using F -tests, using the nominal F -distribution to compute p -values. The viper and AUCell methods do not assign molecules to TFs, but rather calculate enrichment or activity scores that are not counts. We therefore fit Gaussian linear models as implemented in limma using empirical Bayes shrinkage of the variance parameters. This approach was also adopted in Holland et al. [2020] . Here also, the design matrix corresponds to the conditions of interest, and inference is carried out using F -tests. 2.4.3 Gamma-model-based simulation framework In this section, we describe a Gamma-model-based simulation framework. The data-generating mechanism follows the assumptions of our proposed statistical model; transfactor is therefore expected to perform well in this scenario. Specifically, the simulation is based on a Gamma distribution for the TF activities , across all transcription factors in cell type 1, chosen, without loss of generality, as the reference cell type. It requires as input a shape parameter ( h = 1 by default) and a scale parameter ( s = 2 by default) for the Gamma distribution. Given these parameters and a gene regulatory network A, a dataset is simulated as follows. Simulate TF activity for the first (reference) cell type: , independently across all TFs t = 1,…, T . For a subset of TFs chosen to be differentially active across conditions, simulate fold-changes between each cell type c and the reference cell type ( c = 1) from a log-normal ( lN ) distribution: δ tc ∼ lN ( µ = 2.5, σ 2 = .3 2 ), independently across cell types c = 2,…, C . In order to induce both increasing as well as decreasing TF activity between cell types, randomly set sign of log-fold-changes log δ tc using independent Bernoulli(1/2) random variables. For each cell type c = 2,…, C , define . Calculate , under assumptions (A1) and (A2), where and we further assume that α g c = 1. Simulate gene expression measures . We consider a scenario with C = 5 cell types and simulate gene expression measures for 50 cells of each type. The true gene regulatory network consists of G = 500 genes, regulated by T = 80 transcription factors, using only activating edges. For each gene, the number of regulating TFs are sampled from a Binomial distribution with 30 draws and success probability of 0.05. If the Binomial draw equals zero, a new draw is taken until a positive number is obtained. The set of regulating TFs for each gene are then chosen randomly. We simulate a random set of 25% of the TFs to be differentially active between cell types. In order to incorporate prior information on TF activity in transfactor , we let where m g = 100 for all genes by default. This is then used to parameterize the prior Dirichlet distribution. Methods are evaluated using the area under the receiver operating characteristic (ROC) curve, reflecting their capability of recovering the TFs that are differentially active between cell types. All methods are evaluated using the true GRN, as well as GRNs where (a) edges in the true GRN have been removed; (b) edges in the true GRN have been added; (c) ‘noisy’ TFs, i.e., TFs that did not contribute to the data-generating process and regulate a random set of genes, are added. 2.4.4 dyngen simulation framework In a second simulation-based evaluation, we use the dyngen framework to simulate scRNA-seq data. dyngen is a stochastic simulation algorithm for dynamic biological processes, which simulates (nucleic acid and protein) molecular abundances at single-cell resolution and explicitly models key biological processes such as transcription, splicing, translation, and regulation. The simulation framework first constructs a gene regulatory network, including both activating as well as repressing interactions between TFs and target genes, from which molecular abundances are simulated. We can therefore use this framework to evaluate TF activity estimation, using the simulated mRNA abundances. As in the Gamma-model-based simulation framework, we benchmark several settings by adding noisy edges, removing true edges, and adding noisy TFs to the true GRN, while checking how well methods can recover the data-generating TFs. The data-generating GRN was constructed from a linear backbone in dyngen and consists of ∼ 600 genes regulated by 100 to 200 TFs; exact numbers vary between simulated datasets. As prior information, we use the ‘strength’ parameter defined in dyngen , which is a positive number used in the data-generating scheme that approximates the strength of the regulatory interaction between a TF and its target genes. Note, however, that this prior information is only a proxy for the true strength of the edges, as several other parameters are at play in this more complex simulation scenario, and the prior information may therefore be considered more noisy, as compared to the Gamma-model-based simulation framework. 2.4.5 Yeast knock-out case study The yeast dataset was previously published in Jackson et al. [2020] and is provided in the Supplementary Material accompanying that paper. The dataset consists of 6, 829 genes and 38, 225 cells, representing 12 gene deletion mutants across 11 different environmental conditions; see Results for more details. The gene regulatory networks estimated as part of the original publication are also provided in the Supplementary Material. We use the global network aggregated across all environmental conditions (i.e., signed network.tsv file in the Supplementary Material of the original paper) to estimate TF activity. The network consists of 2, 445 genes and 129 TFs, with 1, 052 repressing and 11, 176 activating edges. TF activity estimation was performed using the Poisson model allowing for repressions, but not enforcing sparsity. The heatmap of TF activity is created by standardizing estimated TF activities for each TF separately to have a mean of zero and a variance of one across cells. The dendrograms are created using the Ward clustering method [ Ward, 1963 ] based on Euclidean distances, and the heatmap is visualized using the pheatmap R package. Clustering is performed by ‘cutting’ the dendrogram at a certain height in order to obtain the prespecified number of clusters. For the diauxic shift analysis, we verify the activity of TFs identified in Murphy et al. [2015] as being ‘late-induced’, ‘early-induced’, or ‘early-repressed’ according tof analysis, we use the review paper by the TF enrichment analysis performed by the original authors. For the nitrogen metabolism analysis, we use the review paper by Conrad et al. [2014] to identify four TFs that should be upregulated during a lack of preferred nitrogen sources. To investigate the stability of the results with respect to the cells being sequenced, we perform a bootstrap analysis, where we randomly resample cells with replacement to reach an equally large dataset, but where the number of cells for each environmental condition or genotype may vary across bootstrap samples. TF activities are then estimated on each of 30 bootstrapped dataset using the same procedure as on the full dataset. To incorporate confidence for each edge in the GRN, we use the inferelator method that was also used in [ Jackson et al.2020 ]. While the global network in Jackson et al. [2020] was derived using the AMuSR implementation, here we estimate the networks using the Bayesian best subset regression (BBSR) approach implemented in inferelator as it was more readily available for our purpose. We use inferelator’s built-in implementation to bootstrap cells and re-estimate the GRN five times. The accompanying confidence scores for each edge present in the original GRN are leveraged as prior information when estimating TF activity using the hierarchical Poisson model with known Dirichlet prior, where the alphaScale parameter is set to 1, i.e., S = 1 following the notation in the paragraph on ‘Scaling γ gt ’ in Supplementary Methods. This means that we are giving equal weight to the information provided by the prior as well as the data. Finally, we want to compare transfactor to viper and AUCell on equal grounds, as transfactor may be considered to have an a priori advantage by using the groupings of cells via the design matrix, which is not used for viper and AUCell since they estimate TF activity per cell separately. To avoid a possible competititve advantage, we estimate TF activity by grouping sets of at most 10 cells within each environmental condition. For transfactor , we use these cell sets to specify the design matrix. For viper and AUCell , we sum the counts for each cell within each set, and use these sums to estimate TF activity. Such an analysis allows an unbiased comparison between transfactor , which is supervised via the design matrix, and AUCell and viper , which are unsupervised in that sense. 2.4.6 Olfactory epithelium case study The olfactory epithelium dataset was previously published in Brann et al. [2020] . Here, we only focus on the cells that were assigned to the neuronal lineage. Cells belonging to this lineage were selected using the cell assignment procedure in tradeSeq, based on a trajectory fitted using Slingshot [ Van den Berge et al., 2020 , Street et al., 2018 ]. After subsetting, the dataset consists of 14, 618 genes and 6, 810 cells. The gene regulatory network was estimated using grnboost2 , implemented in SCENIC 0.10.2 [ Aibar et al., 2017 ]. The estimated GRN contains 7, 863 genes, 262 transcription factors, and 25, 896 edges. The median number of genes regulated by a TF is 27. TF activity is estimated using the Poisson model with sparse initialization. We partition the cells into 20 bins, each containing ∼ 340 cells, according to the cells’ pseudotime inferred by Slingshot . We use this binning factor to construct a design matrix for estimating TF activity with transfactor . The contribution of a TF to a cell’s overall gene expression is calculated (using the formula in Equation (21) in Supplementary Methods) and scaled by the cell’s total gene expression measure, i.e., For gene set enrichment analysis on the three TF groupings as obtained in Results, we select, for each TF, the genes for which that TF is found to produce most of the molecules as compared to the other TFs regulating that gene (i.e., the genes g for which the TF has the highest π gtc among all t ). This set of genes for all TFs in a group is then used for gene set enrichment analysis using gProfiler [ Raudvere et al., 2019 ], with Ensembl annotation v102, as can be accessed using https://biit.cs.ut.ee/gprofiler_archive3/e102_eg49_p15/gost . For dimensionality reduction, we perform principal component analysis (PCA) using scater [ McCarthy et al., 2017 ] to derive the top 10 principal components on the output of each method. Since our method outputs counts, we log-transform them before PCA. The top 10 PCs are then used as input for uniform manifold approximation and projection (UMAP) dimensionality reduction [ McInnes et al., 2018 ] with parameter min dist=0.8. 3. Results 3.1 Simulation studies 3.1.1 Gamma-model-based simulation Estimation of TF activity parameters In the gamma-model-based simulation framework, we know the true TF activities and their corresponding contributions to gene expression. Hence, we can evaluate how close the estimates from transfactor are to the true parameter values using the mean squared error (MSE) (Supplementary Feigure 1). Note that this is not possible for viper and AUCell , as their output is on a different scale and therefore non-interpretable in terms of TF activity in the way we simulate it. The poisson and dirMultEBayes models perform similarly in terms of MSE across all conditions, while providing prior information using the dirMult model results in more accurate estimation, especially when the GRN contains many false edges and truly non-active TFs. For all models, the sparse initialization has no influence on the MSE, likely because the truly non-active TFs have very low either way; therefore, setting them to zero has negligible impact on the MSE, which is dominated by truly active TFs with high . Detection of differential TF activity Next, we turn to assessing the performance of TF activity inference methods in terms of their ability to detect differentially active TFs. We use the area under the receiver operating characteristic curve (AUC) as a performance measure; a high AUC corresponds to a high sensitivity and specificity for identifying the TFs truly involved in the system under study. We first discuss the scenario where we know and use the true data-generating GRN (middle panel in top row of Figure 1 ). In this scenario, using a sparse initialization of the EM algorithm provides no benefits, as expected. While all methods perform generally well, the poisson model performs best, as does AUCell . Note that, while it may seem counterintuitive at first that the dirMult model performs somewhat worse, this can be explained by the fact that the prior distribution is shared across cell types, effectively leading to moderation of the TF activity estimates by the prior distribution across cell types. Since the performance measure used here relates to differential TF activity, this moderation can lead to an attenuation of differences in TF activity, hence loss of power, whenever the true TF activities are different across cell types. Download figure Open in new tab Figure 1: Gamma-model-based simulation study results. Eight methods were evaluated on each of six simulated datasets. For each method and dataset, the area under the receiver operating characteristic curve for detecting the TFs that are differentially active between cell types was computed. For each method, the six AUCs were summarized using a boxplot. The different panels represent different GRNs used as input to each method. The first row corresponds to a GRN where only TFs that were involved in the data-generating process are used. For rows 2–4, respectively, 80, 160, and 320 ‘irrelevant’ TFs (i.e., TFs that do not contribute to the data-generating process and regulate a random set of genes) were added to the GRN. The different columns represent removing or adding edges to the true GRN. In the middle column, we use the true GRN. Moving to the left (right) of the middle column, we remove (add) edges to the true GRN. Panel subtitles indicate the number of irrelevant TFs added to the GRN and edges added/removed from the GRN. Next, we evaluate a more realistic scenario, where we do not know the true GRN (see Methods for details). First, we consider the scenario where edges are either removed from or added to the true GRN (i.e., the first row in Figure 1 ). As edges are removed, all methods drop in performance, but both the poisson model and AUCell still retain good performance. viper is the only method for which performance declines severely. As more false positive edges are added to the GRN, there is a markedly larger drop in performance for most methods. As could be expected, the dirMult methods are capable of using the prior information to prioritize true edges in the GRN and thus perform best, while the dirMultEBayes methods perform worst. Among the competing methods, viper seems more robust to false positive edges than AUCell . We also consider the scenario where ‘irrelevant’ TFs are added to the GRN. In general, we see that sparse initialization can become beneficial here as, indeed, the irrelevant TFs are identified and their activities are shrunken to zero. This is especially apparent in the setting where the GRN has both many irrelevant TFs and many false positive edges, i.e., the bottom right panels of Figure 1 . Note that this is likely also the setting that corresponds most closely to a real data analysis, when using a GRN that could, for example, be downloaded from a database and not necessarily be appropriate for the specific context. The dirMult methods perform consistently very well, since they are able to use prior information to prioritize true edges and data-generating TFs. In settings where the GRN can be expected to contain a lot of noise, but there is no prior information available, the poisson method with sparse initialization performs well in general. For the differential activity (DA) analysis based on estimated TF activities from transfactor , we have made use of count models as implemented in glmGamPoi [ Ahlmann-Eltze and Huber, 2020 ], while we have used limma [ Smyth, 2004 ] for viper and AUCell , as their output cannot be considered to be count data. In order to prevent such differences in post-processing from confounding the comparison of transfactor , AUCell , and viper , we also provide a version of Figure where limma is used to perform differential activity analysis for all methods (Supplementary Feigure 2) and the voom transformation approach is applied for all transfactor methods [ Law et al., 2014 ]. While in general the performance of transfactor is lower, most of the results are qualitatively similar. The dirMult model, however, seems to specifically suffer from the change in DA analysis methodology. Finally, AUCell performs better than transfactor when the GRN is missing true edges. Taken together, these results suggest that, if reliable prior information is available, it is beneficial, as expected, to use the dirMult model. However, in real data settings, we generally don’t have as accurate prior information as we have here. When the GRN can be expected to contain both false positive edges and TFs that are not involved in the process being studied, the sparse initialization step helps, and the poisson model could be a good choice. In cases where we expect the GRN to be missing true edges, AUCell also performs well. 3.1.2 dyngen simulation We also evaluate TF activity estimation using the dyngen simulation framework [ Cannoodt et al., 2020 ]. We simulate six datasets, each of which can be represented by a trajectory consisting of a single lineage in a low-dimensional embedding (see Supplementary Feigure 3a, for an example). Note that, since dyngen uses a gene regulatory network to simulate the data (Supplementary Feigure 3c–d), we know the true gene regulatory network that is driving the development represented by this trajectory. Each dataset consists of ∼ 700 cells, for each of which we are observing 100–200 TFs, ∼ 500 target genes, and ∼ 100 housekeeping genes (that are not regulated by any TFs). The GRN used to generate the data includes activating as well as repressing edges, additionally allowing us to benchmark our approach to account for repressions in the GRN. The prior information in this setting is more noisy as compared to the Gamma-model-based simulation framework, due to a more complex simulation scenario. First, evaluating all methods using the true GRN as input (middle panel in top row of Figure 2 ), it is clear that transfactor outperforms competing methods viper and AUCell . As for the gamma-model-based simulation of Section 3.1.1 , sparse initialization has no apparent impact on the performance of the transfactor methods if the TF activities are in fact not sparse. Accounting for repressions in addition to the sparse initialization does not seem to increase performance in this setting. When true edges are removed from the GRN (left panels in top row of Figure 2 ), all methods similarly drop in performance, while AUCell remains rather stable. As false edges are added to the GRN (right panels in top row of Figure 2 ), using the available prior information on the strength of TF-gene interactions from the dyngen simulation framework provides stability to the results, and accounting for repressions can become more beneficial. Overall, our transfactor methods perform better than AUCell and viper . Download figure Open in new tab Figure 2: dyngen simulation study results. Eleven methods were evaluated on each of six simulated datasets. For each method and dataset, the area under the receiver operating characteristic curve for detecting the TFs that are differentially active between the cell types was computed. For each method, the six AUCs were summarized using a boxplot. The different panels represent different GRNs used as input to each method. The first row corresponds to a GRN where only TFs that were involved in the data-generating process are used. For rows 2–4, respectively, 131, 262, and 524 ‘irrelevant’ TFs (i.e., TFs that do not contribute to the data-generating process and regulate a random set of genes) were added to the GRN. The different columns represent removing or adding edges to the true GRN. In the middle column, we use the true GRN. Moving to the left (right) of the middle column, we remove (add) edges to the true GRN. Panel subtitles indicate the number of irrelevant TFs added to the GRN and edges added/removed from the GRN. As more irrelevant TFs are added to the GRN (rows 2–4 in Figure 2 ), sparse initialization results in larger gains, especially when many false positive edges are present in the GRN. Note that the sparse initialization has no real impact on the dirMult model, since it is overridden by the prior information. As edges are removed from the true GRN as well as a large number of irrelevant TFs are added, the poisson , dirMult, and dirMultEBayes methods start performing similarly. Note that when many false positive edges are added to the GRN as well as many irrelevant TFs, the dirMult method is inferior to the poisson method with sparse initialization. This is due to multiple reasons. First, the prior includes information on edges from irrelevant TFs in the GRN, which overrides the sparse initialization (see Methods), hence losing its performance benefit. It is thus more appropriate to compare the dirMult methods to the poisson method without sparse initialization. Second, the prior parameters for noisy and true edges are often of similar magnitude, due to the strength parameters provided by dyngen . The effect observed is therefore partially the result of a less discriminative prior; see Supplementary Feigure 4 for the results when using a prior that is always higher for true vs. false edges. Third, the prior parameters here are only a proxy for their true value, while in the Gamma-model-based simulation study their true value was known. We again investigate performance for all methods when exclusively using limma (-voom) for differential activity analysis. In this simulation study, we notice little influence of the downstream DA approach on methods’ performance (Supplementary Feigure 5). In summary, these results largely confirm those obtained with the Gamma-model-based simulation framework, but additionally show that accounting for true repressing interactions can be beneficial. If reliable prior information on the interaction strength is available, it would be a good idea to use it, however, noisy prior information can also be detrimental, especially if it overrides the sparse initialization step that can be quite powerful. If such information is not available, methods based on the Poisson and empirical Bayes hierarchical Poisson models are powerful alternatives. 3.2 Case studies 3.2.1 Yeast knock-out dataset We analyze a scRNA-seq dataset on budding yeast ( Saccharomyces cerevisiae ) [ Jackson et al., 2020 ]. The dataset encompasses 12 transcriptionally-barcoded gene deletion mutants, each investigated in 11 different environmental conditions, comprising a total 38, 225 single cells. As noted by the authors, the different environmental conditions represent the largest source of variability in the dataset, with the effects of mutations hardly noticeable in reduced-dimensional space [ Jackson et al., 2020 ]. The environmental conditions are summarized in Table 1 . The original publication used the diversity of transcriptional states, in combination with known regulations, to construct a global gene regulatory network across mutants and environmental conditions, consisting of 2, 445 genes and 129 TFs, including 1, 052 repressing and 11, 176 activating edges. We use this GRN to estimate TF activity, accounting for repressions, using our Poisson model. The design matrix accommodates a possible interaction of genotype (gene deletions) and environmental condition. Given the careful construction of the GRN using a range of transcriptional states and prior information relevant to the dataset, we do not enforce sparsity. View this table: View inline View popup Download powerpoint Table 1: Yeast knock-out case study. Environmental conditions. The heatmap in Figure 3a displays the estimated mean TF activities for each mutant in each environmental condition for each of the 129 TFs from the GRN, where, for each TF, activities are standardized to have zero mean and unit variance across cells. As in the original publication, samples cluster mainly by environmental condition. The estimated TF activities are useful summaries of the data, as indicated by their effectiveness at recovering the environmental conditions, here, by hierarchical clustering followed by cutting the dendrogram at a certain height such that the number of clusters corresponds to the number of environmental conditions. The adjusted Rand index (ARI) of such a TF-activity-based clustering corresponds to 0.83, as compared to an ARI of only 0.53 when clustering samples according to the average expression of those same TFs (here, the ARI measures the similarity between the clusterings of samples based on TF activity or expression measures and the reference grouping of samples based on environmental condition). The heatmap of the standardized estimated TF activities highlights some remarkable patterns, in particular, the condition-specific higher activity of TFs and transcriptional states for some conditions that seem to be governed by a limited number of most active TFs, while other TFs show relatively moderate to low activities, across most mutants. Examples are the YPD and CStarve conditions, amongt others. In the YPD condition, where all nutrients are available, a group of seven TFs ( IFH1, SFP1, FKH2, HMS1, HCM1, URC2 , and MIG1 ) cluster together as being most active in that condition (purple box in Figure 3a ). This group of TFs includes IFH1 , a TF known to regulate ribosomal protein (RP) genes and recruited in optimal growth conditions [ Schawalder et al., 2004 ]. SFP1 is known to be essential for normal yeast cell growth [Blumberg and 1991], and also regulates RP genes [ Jorgensen et al., 2004 ], while also contributing to G2/M transitions in the mitotic cell cycle [ Xu and Norris, 1998 ]. FKH2 is similarly involved in cell cycle dynamics, being a known activator of replication origins, and regulates expression of G2/M phase genes. Not much is known about HMS1 , however, overexpression of the TF results in hyperfilamentous growth, hence it could be required for yeast cell growth [ Lorenz and Heitman, 1998 ]. HCM1 drives S-phase activation of genes involved in chromosome segregation [ Pramila et al., 2006 ]. is involved in glucose repression, whereby cells grown on glucose repress the expression of a large number of genes that are required for the metabolism of alternative carbon sources [ Carlson, 1999 ]. Finally, there does not seem to be much known about URC2 , a TF that is possibly involved in uracil catabolism [ Andersen et al., 2008 ], a function that would align with those of other TFs in this group. Download figure Open in new tab Figure 3: Yeast knock-out case study. (a) Heatmap of estimated TF activities from transfactor , where rows correspond to TFs and columns to cell groups based on condition and genotype. The pseudocolor correspond to TF activity, standardized to have zero mean and unit variance across cell groups. (b) TF activity estimates from transfactor for nine TFs involved in the diauxic shift. Each color represents a different gene deletion mutant. (c) TF activity estimates from transfactor for four TFs involved in nitrogen metabolism. Each color represents a different gene deletion mutant. (d) UMAP dimensionality reduction plots for the gene expression counts from aggregated cells (‘Summed counts’, see Methods) and corresponding TF activity estimates using viper , AUCell , and transfactor . Cells are colored according to environmental condition. Aside from the initial exploratory data analysis, we also verify if transfactor can capture regulators known in the literature to be associated with changes caused by particular environmental conditions. Murphy et al. [2015] study protein abundances during the diauxic shift using a proteomics time series experiment. They identify protein lists that are either induced or repressed, early or late in the diauxic shift, and subsequently use the YEASTRACT transcription factor enrichment tool [ Teixeira et al., 2023 ] to assess which transcription factors might drive the observed differences in protein abundance. This procedure identifies a total of nine TFs, seven of which should be upregulated in the diauxic shift and two which should be downregulated. Our estimated TF activities confirm the findings of the proteomic study, and show very good agreement of up- and downregulation of TFs in the diauxic shift, as compared to the control (YPD) condition ( Figure 3b ). One notable exception is the BAS1 TF, where the up- or downregulation is variable across different mutants. For HAP3 , we find very low TF activity in both conditions. The dataset also includes four environmental conditions that are related to nitrogen limitation, where the preferred nitrogen source (ammonium sulfate) is replaced by alternative sources. Conrad et al. [2014] review nutrient sensing and signaling in yeast. They discuss that, in the absence of a preferred nitrogen source, four TFs ( RTG1, RTG3, GLN3 , and GAT1 ) are activated, while in the presence of a preferred nitrogen source these are deactivated through (hyper)phosphorylation. We therefore compare TF activity of all four nitrogen-limitation conditions with the YPD condition, and indeed detect increased TF activity in all nitrogen conditions for most of these TFs, confirming what is known in the literature ( Figure 3c ). While confirming known biology is important, the results should also be stable with respect to the cells being sequenced. We therefore performed a bootstrap analysis, by randomly resampling cells with replacement to get an equally large dataset, and re-estimating TF activity to see if our conclusions still stand. The changes in TF activity for both the diauxic shift as well as the nitrogen limitation TFs are stable across bootstrap samples (Supplementary Feigures 6–7). Furthermore, we assess the stability of the results when incorporating information on the confidence of the edges in the GRN. To this end, we re-estimate the GRN across five different bootstrap samples using the inferelator method (see Methods), which the authors of the study used to infer their GRN, and leverage the confidence scores (i.e., the probability each edge is present across the bootstrap samples) for each edge as prior information in the TF activity estimation using the hierarchical Poisson model. Incorporating the uncertainty on GRN edges results in TF activities that equally well cluster the samples according to environmental condition, with an ARI of 0.82. The biological findings on the diauxic shift and nitrogen metabolism are also confirmed, with the exception that the BAS1 TF is now less active in the diauxic condition as compared to the control (Supplementary Feigure 8). We also compare transfactor with viper and AUCell in terms of recovering the biological structure (i.e., the different environmental conditions) of the data upon TF activity estimation. Both viper and AUCell are unsupervised methods, assessing TF activity for each cell separately, while transfactor estimates TF activity separately within each of a predefined set of cell types. To be able to compare all methods on an equal footing, we first group cells in random sets of 10 within each environmental condition, and use this grouping as design matrix for transfactor , while for viper and AUCell we add counts for these sets of cells and then assess TF activity estimation. This analysis strategy ensures that all methods use the same information and avoids a priori competitive advantages. Figure 3d shows UMAP representations of the counts and corresponding TF activity estimates from all methods. While viper and transfactor recover much of the biological variability, AUCell fails at clearly separating many environmental conditions. This is also reflected in the ARI, which is 0.32 for both viper and transfactor , and 0.21 for AUCell . Finally, one may wonder whether the results obtained in this case study may have similarly been obtained without estimating TF activity, and solely considering TF expression instead. We have already seen that TF activity estimates better reflect the biological structure of the dataset, as noted by the increased ARI when comparing activity-based (vs. expression-based) clusters to environmental conditions. The confirmation of increased and decreased TF activity in the diauxic shift condition are much less clear when looking at TF expression rather than TF activity. Indeed, the expression of BAS1 and MSN2 conflicts with what is known in the literature. Further, the two TFs that should decrease in activity ( FHL1 and RAP1 ) now no longer consistently do so in terms of expression (Supplementary Feigure 9a). Similarly, results on the nitrogen limitation condition show less agreement with the literature when looking at TF expression as compared to TF activity (Supplementary Feigure 9b). In summary, this case study shows that TF activity estimates are relevant and useful, allowing novel insight in terms of TF regulation. It additionally demonstrates the flexibility of our method to handle different situations with varying levels of information on the GRN. 3.2.2 Olfactory epithelium dataset The olfactory epithelium (OE) is responsible for the perception of smell and sending these molecular signals to the brain. It also serves as a barrier between the outer world and the brain, requiring neuronal cells to be continuously replaced through a source population of progenitor cells. In Brann et al. [2020] , the OE is chemically injured to trigger regeneration and repair of the tissue through activation of stem cells known as horizontal basal cells (HBC). To develop into neuronal cells, HBCs first differentiate into globose basal cells (GBC), which mitotically divide toensure continuous neuronal regeneration [ Fletcher et al.,, 2017 ]. These GBCs in turn develop into immature neurons and finally mature neurons [ Fletcher et al., 2017 , Gadye et al., 2017 ]. The trajectory is shown in Figure 4a . Download figure Open in new tab Figure 4: Olfactory epithelium case study. (a) UMAP dimensionality reduction of the gene expression data, colored according to cell type. A clear differentiation path from activated horizontal basal cells (HBC*) to neuronal cells is observed. (b) Heatmap of estimated TF activities from transfactor for the 235 TFs found to be involved in the neuronal lineage. Each row represents a TF and the different columns represent 20 cell groups of approximately equal size, binned according to their pseudotime values. The pseudocolor correspond to TF activity, standardized to have zero mean and unit variance across bins. Three clusters are identified based on their activity patterns. (c) Activity patterns of the top 4 transcription factors, where the top list was derived using the scaled Euclidean distance measure (see Equation (4) in Methods). The y-axis represents the contribution of a TF to the entire observed gene expression in each cell ( Equation (19) ), e.g., if the contribution equals 0.1, then we estimate that that TF was involved in the production of 10% of all molecules present in the cell. The x-axis represents the 20 pseudotime bins. Note that panels are on different y-axis scales. (d) UMAP reduced-dimensionality plots based on estimated TF activity for four methods. Cells are colored according to their pseudotime bin, from dark blue (first bin) to red (last bin). Note that panels are on different scales. co-regulated modules. Alternatively, community detection algorithms (e.g., [ Li et al., 2020 ]) could be used. TF activity is estimated using the transfactor poisson model with sparse initialization. Out of a total of 262 TFs that are included in the GRN estimated by SCENIC [ Aibar et al., 2017 ], our model finds that 27 do not actively contribute to the process, i.e., zero molecules are associated with these TFs. A heatmap of the standardized activities of the contributing TFs shows three distinct groups of TFs which can be discriminated based on their activity patterns ( Figure 4b ). These groups are respectively most active in the HBC, GBC, and neuronal stages, recapitulating the biological knowledge on these three previously-identified cell types [ Fletcher et al., 2017 , Gadye et al., 2017 ]. The relevance of the grouping is confirmed by gene set enrichment analysis on the genes found to be regulated by these TF groups, which identifies significant biological processes involved in tissue development, cell cycle, and cell projection for, respectively, the HBC, GBC, and neuronal groups. The top enriched gene sets for each group are provided in Supplementary Tables 1–3. High-throughput studies typically focus on differential expression analysis to investigate gene regulation. Here, we identify differentially expressed genes along the neuronal lineage using tradeSeq [ Van den Berge et al., 2020 ] and we find 6, 977 DE genes at a 5% nominal false discovery rate (FDR) level; a long list of genes that is cumbersome to interpret. Often, one tries to shorten the list of DE genes by imposing a fold-change cut-off (e.g., McCarthy and Smyth[2009b]). However, testing against a fold-change cut-off of 2, still results in 4, 752 DE genes. In comparison, again using tradeSeq, we find that 209 TFs are differentially active along the neuronal lineage, 135 of which are still differentially active after testing against the same fold-change threshold, leading to a much more parsimonious interpretation of the dataset. We use a scaled Euclidean distance-based ranking of transcription factors (see Methods) to order them according to their contribution to the observed differences in gene expression. Here, we focus on the top 4 transcription factors, which collectively explain 28% of the total scaled Euclidean distance across all TFs. The TF activities, scaled by the the total activity across all TFs, are shown in Figure 4c . The top transcription factor, Hdac2 , is a member of the Hdac family, which is interestingly known to be essential in the response of neurons following injury [ Cho and Cavalli, 2014 ]. Rfx3 is required for cilia formation in neurogenesis [ Choksi et al., 2014 , is known to be involved in the injury, activation, and differentiation of stem cells [ Okada et al., 1996 ,, 2015]. Fos and Weems, 1999], and is indeed most active at the early stages of development in the lineage. Ezh2 is most active at the GBC stage; previous work confirms that Ezh2 protein activity is indeed found to be highest in the GBC stage [ Goldstein et al., 2018 ]. In summary, the literature confirms that our top transcription factors are involved in olfactory neuron maturation and regeneration upon injury. On a larger scale, we check how well our method, but also viper and AUCell , are able to recover the global structure of the data following TF activity estimation. We therefore run viper and AUCell on all cells using the same GRN estimated by SCENIC . Note, however, that in this analysis our method makes explicit use of the grouping of cells into 20 equally-sized bins according to pseudotime, while viper and AUCell are completely unsupervised, and we therefore expect our method to already have a competitive advantage to recover the structure defined by the cell groupings. We perform dimensionality reduction on the output of each method, by using the top 10 principal components as input to UMAP [ McInnes et al., 2018 ]. The global structure is generally recovered in reduced-dimensional space for all methods, where a differentiation path of cells corresponding to the pseudotime of the trajectory can be observed ( Figure 4d ). Finally, we investigate the robustness of our results to an analysis where we do not use the sparse initialization step. All four top TFs discussed previously are still the top four. In general, the rankings, especially the top rankings, are highly correlated. Performing the same dimensionality reduction on the non-sparse output again leads to recovering the global structure of the data ( Figure 4d ). 4 Discussion We have proposed the transfactor method to infer transcription factor activity by deconvolving transcription factor-specific gene expression measures from observed overall gene expression measures, based on a Poisson model for read counts. A hierarchical Bayes version of the model allows the incorporation of prior information on TF activity. The method is implemented in the open-source software package transfactor available at https://github.com/koenvandenberge/transfactor . The transfactor method relies on a gene regulatory network as input, and we have found the GRN to have a large impact on downstream results. However, it is currently unclear how to assess whether a GRN, either obtained from databases or estimated empirically from the data, is appropriate for a specific dataset. For example, one would perhaps be able to check this by examining whether co-regulated genes are indeed correlated in terms of their measured gene expression. Using TF activity instead of gene expression in downstream analyses has proven to be useful in multiple ways. By substantially reducing the dimensionality of the dataset, use of TF activities should increase the signal-to-noise ratio and allow for easier biological interpretation. In our case studies, TF activities have shown to be better than gene expression measures at grouping cells with respect to their biological condition, and the obtained changes in transcription factor activity match what is known from the literature. The estimated TF activities allow one to prioritize TFs based on their estimated contributions to observed gene expression changes. However, a comparison of TF activity across groups of cells or along lineages of a trajectory, as done in this manuscript, effectively ignores the inherent uncertainty on the estimates. Further work should be directed to estimating the uncertainty on the estimated TF activities. For example, in transcript-level quantification of RNA-seq, where an EM algorithm is adopted to derive transcript abundances, the uncertainty is quantified by bootstrapping the data and re-running the EM algorithm [ Bray et al., 2016 , Patro et al., 2017 ]. Other approaches such as described in Louis [1982] allow for uncertainty estimation to be embedded in the framework of the EM algorithm [ Bray et al., 2016 , Patro et al., 2017 ], however, significantly increasing its computational burden. One of the core assumptions of our model is that, if a TF is active, it is so across all of its target genes. This is surely a simplification of reality, where transcription factors may vary in their effective DNA binding according to specific sequence preferences and therefore between different target genes. An interesting direction would be to weigh the edges of the GRN using models that are capable of computationally estimating binding probabilities to specific target sequences (in this case, gene promotors), such as the DeepBind or BPNet models [ Alipanahi et al., 2015 , Avsec et al., 2021]. In addition, integration with other data modalities, such as the single-cell assay for transposase-accessible chromatin (scATAC-seq), could be relevant in order to detect plausible binding regions and use TF-specific binding probabilities for each of these regions. Finally, one may also consider a different parameterization of the statistical model, e.g., no longer the multiplicative model that we propose here, but a more flexible parameterization. The sparse initialization step of the EM algorithm relies on an equal partitioning of molecules to the TFs that could have produced them. This leaves the sparse initialization step vulnerable to errors, possibly leading to both false positives and false negatives. To avoid this, one could first run the EM algorithm for a few iterations (using the same starting values as usual), after which the obtained estimates could be used in the sparse initialization step. The subsequent sparse solution can then be used as input to the following iterations. While in our manuscript we have focused on exploring the estimated TF activities as the main tool for interpreting the datasets, our method can also be used for more detailed analysis of co-regulation. For example, for each cell type, one could construct a weighted network with the same structure as the GRN A, but where the edges have weights equal to . If one considers A as a bipartite graph, recent biclustering techniques (e.g., [ Moran et al., 2021 ]) could be applied to simultaneously group genes and transcription factors together, in order to identify Acknowledgements The authors thank John Ngai for sharing the dataset that was used in the olfactory epithelium case study and for providing feedback on the obtained biological results. References ↵ Greg Finak , Andrew McDavid , Masanao Yajima , Jingyuan Deng , Vivian Gersuk , Alex K. Shalek , Chloe K. Slichter , Hannah W. Miller , M. Juliana McElrath , Martin Prlic , Peter S. Linsley , and Raphael Gottardo . MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data . Genome Biology , 16 ( 1 ): 278 , 12 2015 . ISSN 1474-760X . doi: 10.1186/s13059-015-0844-5 . URL http://genomebiology.com/2015/16/1/278 . OpenUrl CrossRef PubMed ↵ Koen Van den Berge , Hector Roux de Bézieux , Kelly Street , Wouter Saelens , Robrecht Cannoodt , Yvan Saeys , Sandrine Dudoit , and Lieven Clement . Trajectory-based differential expression analysis for single-cell sequencing data . Nature Communications , 11 ( 1 ): 1201 , 12 2020 . ISSN 2041-1723 . doi: 10.1038/s41467-020-14766-3 . URL http://www.nature.com/articles/s41467-020-14766-3 . OpenUrl CrossRef ↵ GK Smyth . Linear models and empirical bayes methods for assessing differential expression in microarray experiments . Statistical applications in genetics and molecular biology , 3 ( 1 ): 1 – 26 , 2004 . URL http://www.degruyter.com/view/j/sagmb.2004.3.1/sagmb.2004.3.1.1027/sagmb.2004.3.1.1027.xml . OpenUrl ↵ José L. McFaline-Figueroa , Andrew J. Hill , Xiaojie Qiu , Dana Jackson , Jay Shendure , and Cole Trapnell . A pooled single-cell genetic screen identifies regulatory checkpoints in the continuum of the epithelial-to-mesenchymal transition . Nature Genetics , 51 ( 9 ): 1389 – 1398 , 9 2019 . ISSN 1061-4036 . doi: 10.1038/s41588-019-0489-5 . URL http://www.nature.com/articles/s41588-019-0489-5 . OpenUrl CrossRef PubMed ↵ Davis J McCarthy and Gordon K Smyth . Testing significance relative to a fold-change threshold is a TREAT . Bioinformatics , 25 ( 6 ): 765 – 771 , 3 2009a . OpenUrl CrossRef PubMed Web of Science ↵ Samuel A. Lambert , Arttu Jolma , Laura F. Campitelli , Pratyush K. Das , Yimeng Yin , Mihai Albu , Xiaoting Chen , Jussi Taipale , Timothy R. Hughes , and Matthew T. Weirauch. The Human Transcription Factors . Cell , 172 ( 4 ): 650 – 665 , 2 2018 . ISSN 10974172 . doi: 10.1016/J.CELL.2018.01.029/ATTACHMENT/ OpenUrl CrossRef PubMed EDE37821-FD6F-41B7-9A0E-9D5410855AE6/MMC2.XLSX . URL http://www.cell.com/article/S0092867418301065/fulltext http://www.cell.com/article/S0092867418301065/abstract https://www.cell.com/cell/abstract/S0092-8674(18)30106-5 . ↵ Adam Frankish , Mark Diekhans , Anne Maud Ferreira , Rory Johnson , Irwin Jungreis , Jane Loveland , Jonathan M. Mudge , Cristina Sisu , James Wright , Joel Armstrong , If Barnes , Andrew Berry , Alexandra Bignell , Silvia Carbonell Sala , Jacqueline Chrast , Fiona Cunningham , Tomás Di Domenico , Sarah Donaldson , Ian T. Fiddes , Carlos García Girón , Jose Manuel Gonzalez , Tiago Grego , Matthew Hardy , Thibaut Hourlier , Toby Hunt , Osagie G. Izuogu , Julien Lagarde , Fergal J. Martin , Laura Martínez , Shamika Mohanan , Paul Muir , Fabio C.P. Navarro , Anne Parker , Baikang Pei , Fernando Pozo , Magali Ruffier , Bianca M. Schmitt , Eloise Stapleton , Marie Marthe Suner , Irina Sycheva , Barbara Uszczynska-Ratajczak , Jinuri Xu , Andrew Yates , Daniel Zerbino , Yan Zhang , Bronwen Aken , Jyoti S. Choudhary , Mark Gerstein , Roderic Guigó , Tim J.P. Hubbard , Manolis Kellis , Benedict Paten , Alexandre Reymond , Michael L. Tress , and Paul Flicek . GENCODE reference annotation for the human and mouse genomes . Nucleic acids research , 47 ( D1 ): D766 – D773 , 1 2019 . ISSN 1362-4962 . doi: 10.1093/NAR/GKY955 . URL https://pubmed.ncbi.nlm.nih.gov/30357393/ . OpenUrl CrossRef PubMed ↵ Caleb Weinreb , Alejo Rodriguez-Fraticelli , Fernando D. Camargo , and Allon M. Klein . Lineage tracing on transcriptional landscapes links state to fate during differentiation . Science , 367 ( 6479 ), 2 2020 . ISSN 0036-8075 . doi: 10.1126/SCIENCE.AAW3381 . URL https://science.sciencemag.org/content/367/6479/eaaw3381 https://science.sciencemag.org/content/367/6479/eaaw3381.abstract . OpenUrl CrossRef ↵ Amy F Chen , Benjamin Parks , Arwa S Kathiria , Benjamin Ober-Reynolds , Jorg Goronzy , and William J Greenleaf . NEAT-seq: Simultaneous profiling of intra-nuclear proteins, chromatin accessibility, and gene expression in single cells . bioRxiv , page 2021.07.29.454078, 7 2021 . doi: 10.1101/2021.07.29.454078 . URL https://www.biorxiv.org/content/10.1101/2021.07.29.454078v1 https://www.biorxiv.org/content/10.1101/2021.07.29.454078v1.abstract . OpenUrl Abstract / FREE Full Text ↵ Ivan Berest , Christian Arnold , Armando Reyes-Palomares , Giovanni Palla , Kasper Dindler Rasmussen , Holly Giles , Peter-Martin Bruch , Wolfgang Huber , Sascha Dietrich , Kristian Helin , and Judith B Zaugg . Quantification of Differential Transcription Factor Activity and Multiomics-Based Classification into Activators and Repressors: diffTF . Cell reports , 29 ( 10 ): 3147 – 3159 , 12 2019 . ISSN 2211-1247 . doi: 10.1016/j.celrep.2019.10.106 . URL http://www.ncbi.nlm.nih.gov/pubmed/31801079 . OpenUrl CrossRef PubMed ↵ Jonathan D. Rubin , Jacob T. Stanley , Rutendo F. Sigauke , Cecilia B. Levandowski , Zachary L. Maas , Jessica Westfall , Dylan J. Taatjes , and Robin D. Dowell . Transcription factor enrichment analysis (TFEA): Quantifying the activity of hundreds of transcription factors from a single experiment . bioRxiv , page 2020.01.25.919738, 2 2020 . doi: 10.1101/2020.01.25.919738 . URL https://www.biorxiv.org/content/10.1101/2020.01.25.919738v3 . OpenUrl Abstract / FREE Full Text ↵ Mariano J Alvarez , Yao Shen , Federico M Giorgi , Alexander Lachmann , B Belinda Ding , B Hilda Ye , and Andrea Califano . Functional characterization of somatic mutations in cancer using network-based inference of protein activity . Nature Genetics , 48 ( 8 ): 838 – 847 , 8 2016 . ISSN 1061-4036 . doi: 10.1038/ng.3593 . URL nature.com/articles/ng.3593. OpenUrl CrossRef PubMed ↵ Sara Aibar , Carmen Bravo González-Blas , Thomas Moerman , Vân Anh Huynh-Thu , Hana Imrichova , Gert Hulselmans , Florian Rambow , Jean-Christophe Marine , Pierre Geurts , Jan Aerts , Joost van den Oord , Zeynep Kalender Atak , Jasper Wouters , and Stein Aerts . SCENIC: single-cell regulatory network inference and clustering . Nature Methods , 14 ( 11 ): 1083 – 1086 , 11 2017 . ISSN 1548-7091 . doi: 10.1038/nmeth.4463 . URL http://www.nature.com/articles/nmeth.4463 . OpenUrl CrossRef PubMed ↵ Christian H. Holland , Jovan Tanevski , Javier Perales-Patón , Jan Gleixner , Manu P. Kumar , Elisabetta Mereu , Brian A. Joughin , Oliver Stegle , Douglas A. Lauffenburger , Holger Heyn , Bence Szalai , and Julio Saez-Rodriguez . Robustness and applicability of transcription factor and pathway analysis tools on single-cell RNA-seq data . Genome Biology , 21 ( 1 ): 36 , 12 2020 . ISSN 1474-760X . doi: 10.1186/s13059-020-1949-z . URL https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-1949-z . OpenUrl CrossRef PubMed ↵ Andrew J Bannister and Tony Kouzarides . Regulation of chromatin by histone modifications . Cell Research , 21 ( 3 ): 381 – 395 , 3 2011 . ISSN 1001-0602 . doi: 10.1038/cr.2011.22 . URL http://www.nature.com/articles/cr201122 . OpenUrl CrossRef PubMed Web of Science ↵ Can Sönmezer , Rozemarijn Kleinendorst , Dilek Imanci , Guido Barzaghi , Laura Villacorta , Dirk Schübeler , Vladimir Benes , Nacho Molina , and Arnaud Regis Krebs . Molecular Co-occupancy Identifies Transcription Factor Binding Cooperativity In Vivo . Molecular cell , 81 ( 2 ): 255 – 267 , 1 2021 . ISSN 1097-4164 . doi: 10.1016/j.molcel.2020.11.015 . URL http://www.ncbi.nlm.nih.gov/pubmed/33290745 . OpenUrl CrossRef ↵ Nicolas L Bray , Harold Pimentel , Páll Melsted , and Lior Pachter . Near-optimal probabilistic RNA-seq quantification . Nature Biotechnology , 34 ( 5 ): 525 – 527 , 4 2016 . ISSN 1087-0156 . doi: 10.1038/nbt.3519 . URL http://www.nature.com/doifinder/10.1038/nbt.3519 . OpenUrl CrossRef PubMed ↵ Rob Patro , Geet Duggal , Michael I Love , Rafael A Irizarry , and Carl Kingsford . Salmon provides fast and bias-aware quantification of transcript expression . Nature Methods , 14 ( 4 ): 417 – 419 , 3 2017 . ISSN 1548-7091 . doi: 10.1038/nmeth.4197 . URL http://www.ncbi.nlm.nih.gov/pubmed/28263959 http://www.nature.com/doifinder/10.1038/nmeth.4197 . OpenUrl CrossRef PubMed ↵ A.P. Dempster , N.M. Laird , and D.B. Rubin . Maximum Likelihood from Incomplete Data via the EM Algorithm on JSTOR . Journal of the Royal Statistical Society. Series B (Methodological) , 39 ( 1 ): 1 – 38 , 1977 . URL http://www.jstor.org/stable/2984875?seq=1#page_scan_tab_contents . OpenUrl CrossRef Web of Science ↵ Jerzy Neyman . Contribution to the theory of the chi-squared test . Proceedings of First Berkeley Symposium , pages 239 – 73 , 1949 . ↵ Robert Tibshirani . Regression Shrinkage and Selection via the Lasso on JSTOR . Journal of the Royal Statistical Society. Series B (Methodological) , 58 ( 1 ): 267 – 288 , 1996 . URL https://www.jstor.org/stable/2346178?seq=1#metadata_info_tab_contents . OpenUrl CrossRef Web of Science ↵ Jerome Friedman , Trevor Hastie , and Robert Tibshirani . Regularization Paths for Generalized Linear Models via Coordinate Descent . Journal of Statistical Software , 33 ( 1 ): 1 – 22 , 2 2010 . ISSN 1548-7660 . doi: 10.18637/jss.v033.i01 . URL http://www.jstatsoft.org/v33/i01/ . OpenUrl CrossRef PubMed ↵ Constantin Ahlmann-Eltze and Wolfgang Huber . glmGamPoi: Fitting Gamma-Poisson Generalized Linear Models on Single Cell Count Data . bioRxiv , page 2020.08.13.249623, 9 2020 . doi: 10.1101/2020.08.13.249623 . URL https://www.biorxiv.org/content/10.1101/2020.08.13.249623v2 . OpenUrl Abstract / FREE Full Text ↵ Christopher A Jackson , Dayanne M Castro , Giuseppe-Antonio Saldi , Richard Bonneau , and David Gresham . Gene regulatory network reconstruction using single-cell RNA sequencing of barcoded genotypes in diverse environments . eLife , 9 , 1 2020 . ISSN 2050-084X . doi: 10.7554/eLife.51254 . URL https://elifesciences.org/articles/51254 . OpenUrl CrossRef ↵ Joe H. Ward . Hierarchical Grouping to Optimize an Objective Function . Journal of the American Statistical Association , 58 ( 301 ): 236 – 244 , 3 1963 . ISSN 0162-1459 . doi: 10.1080/01621459.1963.10500845 . URL http://www.tandfonline.com/doi/abs/10.1080/01621459.1963.10500845 . OpenUrl CrossRef PubMed J Patrick Murphy , Ekaterina Stepanova , Robert A Everley , Joao A Paulo , and Steven P Gygi . Comprehensive Temporal Protein Dynamics during the Diauxic Shift in Saccharomyces cerevisiae . Molecular & cellular proteomics : MCP , 14 ( 9 ): 2454 – 65 , 9 2015 . ISSN 1535-9484 . doi: 10.1074/mcp.M114.045849 . URL http://www.ncbi.nlm.nih.gov/pubmed/26077900 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC4563728 . OpenUrl Abstract / FREE Full Text ↵ Michaela Conrad , Joep Schothorst , Harish Nag Kankipati , Griet Van Zeebroeck , Marta Rubio-Texeira , and Johan M. Thevelein . Nutrient sensing and signaling in the yeast Saccharomyces cerevisiae . FEMS Microbiology Reviews , 38 ( 2 ): 254 – 299 , 3 2014 . ISSN 1574-6976 . doi: 10.1111/1574-6976.12065 . URL https://academic.oup.com/femsre/article-lookup/doi/10.1111/1574-6976.12065 . OpenUrl CrossRef PubMed Web of Science ↵ David H Brann , Tatsuya Tsukahara , Caleb Weinreb , Marcela Lipovsek , Koen Van den Berge , Boying Gong , Rebecca Chance , Iain C Macaulay , Hsin-Jung Chou , Russell B Fletcher , Diya Das , Kelly Street , Hector Roux de Bezieux , Yoon-Gi Choi , Davide Risso , Sandrine Dudoit , Elizabeth Purdom , Jonathan Mill , Ralph Abi Hachem , Hiroaki Matsunami , Darren W Logan , Bradley J Goldstein , Matthew S Grubb , John Ngai , and Sandeep Robert Datta . Non-neuronal expression of SARS-CoV-2 entry genes in the olfactory system suggests mechanisms underlying COVID-19-associated anosmia . Science advances , 6 ( 31 ): eabc5801 , 7 2020 . ISSN 2375-2548 . doi: 10.1126/sciadv.abc5801 . URL http://www.ncbi.nlm.nih.gov/pubmed/32937591 . OpenUrl FREE Full Text ↵ Kelly Street , Davide Risso , Russell B. Fletcher , Diya Das , John Ngai , Nir Yosef , Elizabeth Purdom , and Sandrine Dudoit . Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics . BMC Genomics , 19 ( 1 ): 477 , 12 2018 . ISSN 1471-2164 . doi: 10.1186/s12864-018-4772-0 . URL https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-018-4772-0 . OpenUrl CrossRef PubMed ↵ Uku Raudvere , Liis Kolberg , Ivan Kuzmin , Tambet Arak , Priit Adler , Hedi Peterson , and Jaak Vilo . g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update) . Nucleic Acids Research , 47 ( W1 ): W191 – W198 , 7 2019 . ISSN 0305-1048 . doi: 10.1093/nar/gkz369 . URL https://academic.oup.com/nar/article/47/W1/W191/5486750 . OpenUrl CrossRef PubMed ↵ Davis J. McCarthy , Kieran R. Campbell , Aaron T. L. Lun , and Quin F. Wills . Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R . Bioinformatics , 33 ( 8 ): btw777 , 1 2017 . ISSN 1367-4803 . doi: 10.1093/bioinformatics/btw777 . URL http://www.ncbi.nlm.nih.gov/pubmed/28088763 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC5408845 https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btw777 . OpenUrl CrossRef PubMed ↵ Leland McInnes , John Healy , and James Melville . UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction . Preprint at https://arxiv.org/abs/1802.03426 , x2 2018 . ↵ Charity W Law , Yunshun Chen , Wei Shi , and Gordon K Smyth . voom: Precision weights unlock linear model analysis tools for RNA-seq read counts . Genome biology , 15 ( 2 ): R29 , 1 2014 . ISSN 1465-6914 . doi: 10.1186/gb-2014-15-2-r29 . URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=4053721&tool=pmcentrez&rendertype=abstract . OpenUrl CrossRef PubMed ↵ Robrecht Cannoodt , Wouter Saelens , Louise Deconinck , and Yvan Saeys . dyngen: a multi-modal simulator for spearheading new single-cell omics analyses . bioRxiv, page 2020.02.06.936971 , 9 2020 . doi: 10.1101/2020.02.06.936971 . URL https://www.biorxiv.org/content/10.1101/2020.02.06.936971v3 . OpenUrl Abstract / FREE Full Text ↵ Stephan B. Schawalder , Mehdi Kabani , Isabelle Howald , Urmila Choudhury , Michel Werner , and David Shore . Growth-regulated recruitment of the essential yeast ribosomal protein gene activator Ifh1 . Nature , 432 ( 7020 ): 1058 – 1061 , 12 2004 . ISSN 0028-0836 . doi: 10.1038/nature03200 . URL http://www.nature.com/articles/nature03200 . OpenUrl CrossRef PubMed Web of Science Hal Blumberg and Pamela Silver . A split zinc-finger protein is required for normal yeast growth . Gene , 107 ( 1 ): 101 – 110 , 10 1991 . ISSN 0378-1119 . doi: 10.1016/0378-1119(91)90302-R . URL www.sciencedirect.com/science/article/abs/pii/037811199190302R . OpenUrl CrossRef PubMed Web of Science ↵ Paul Jorgensen , Ivan Rupes , Jeffrey R Sharom , Lisa Schneper , James R Broach , and Mike Tyers . A dynamic transcriptional network communicates growth potential to ribosome synthesis and critical cell size . Genes & development , 18 ( 20 ): 2491 – 505 , 10 2004 . ISSN 0890-9369 . doi: 10.1101/gad.1228804 . URL http://www.ncbi.nlm.nih.gov/pubmed/15466158 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC529537 . OpenUrl Abstract / FREE Full Text ↵ Z Xu and D Norris . The SFP1 gene product of Saccharomyces cerevisiae regulates G2/M transitions during the mitotic cell cycle and DNA-damage response . Genetics , 150 ( 4 ): 1419 – 28 , 12 1998 . ISSN 0016-6731 . URL http://www.ncbi.nlm.nih.gov/pubmed/9832520 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC1460418 . OpenUrl Abstract / FREE Full Text ↵ M C Lorenz and J Heitman . Regulators of pseudohyphal differentiation in Saccharomyces cerevisiae identified through multicopy suppressor analysis in ammonium permease mutant strains . Genetics , 150 ( 4 ): 1443 – 57 , 12 1998 . ISSN 0016-6731 . URL http://www.ncbi.nlm.nih.gov/pubmed/9832522 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC1460428 . OpenUrl Abstract / FREE Full Text ↵ Tata Pramila , Wei Wu , Shawna Miles , William Stafford Noble , and Linda L Breeden . The Forkhead transcription factor Hcm1 regulates chromosome segregation genes and fills the S-phase gap in the transcriptional circuitry of the cell cycle . Genes & development , 20 ( 16 ): 2266 – 78 , 8 2006 . ISSN 0890-9369 . doi: 10.1101/gad.1450606 . URL http://www.ncbi.nlm.nih.gov/pubmed/16912276 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC1553209 . OpenUrl Abstract / FREE Full Text ↵ Marian Carlson . Glucose repression in yeast . Current Opinion in Microbiology , 2 ( 2 ): 202 – 207 , 4 1999 . ISSN 1369-5274 . doi: 10.1016/S1369-5274(99)80035-6 . URL https://www.sciencedirect.com/science/article/pii/S1369527499800356 . OpenUrl CrossRef PubMed Web of Science ↵ Gorm Andersen , Olof Bjornberg , Silvia Polakova , Yuriy Pynyaha , Anna Rasmussen , Kasper Moller , Anders Hofer , Thomas Moritz , Michael Paolo Bastner Sandrini , Anna-Maria Merico , Concetta Compagno , Hans-Erik Akerlund , Zoran Gojkovic , and Jure Piskur . A Second Pathway to Degrade Pyrimidine Nucleic Acid Precursors in Eukaryotes . Journal of Molecular Biology , 380 ( 4 ): 656 – 666 , 7 2008 . ISSN 0022-2836 . doi: 10.1016/J.JMB.2008.05.029 . URL https://www.sciencedirect.com/science/article/pii/S0022283608005846 . OpenUrl CrossRef PubMed Web of Science ↵ Miguel Cacho Teixeira , Romeu Viana , Margarida Palma , Jorge Oliveira , Mónica Galocha , Marta Neves Mota , Diogo Couceiro , Maria Galhardas Pereira , Miguel Antunes , Inês V Costa , Pedro Pais , Carolina Parada , Claudine Chaouiya , Isabel Sá-Correia , and Pedro Tiago Monteiro . YEASTRACT+: a portal for the exploitation of global transcription regulation and metabolic model data in yeast biotechnology and pathogenesis . Nucleic Acids Research , 51 ( D1 ): D785 – D791 , January 2023 . ISSN 0305-1048 . doi: 10.1093/nar/gkac1041 . URL https://www.sciencedirect.com/science/article/pii/S0022283608005846 . OpenUrl CrossRef PubMed ↵ Russell B. Fletcher , Diya Das , Levi Gadye , Kelly N. Street , Ariane Baudhuin , Allon Wagner , Michael B. Cole , Quetzal Flores , Yoon Gi Choi , Nir Yosef , Elizabeth Purdom , Sandrine Dudoit , Davide Risso , and John Ngai . Deconstructing Olfactory Stem Cell Trajectories at Single-Cell Resolution . Cell Stem Cell , 20 ( 6 ): 817 – 830 , 6 2017 . ISSN 19345909 . doi: 10.1016/j.stem.2017.04.003 . URL http://www.ncbi.nlm.nih.gov/pubmed/28506465 . OpenUrl CrossRef PubMed ↵ Levi Gadye , Diya Das , Michael A Sanchez , Kelly Street , Ariane Baudhuin , Allon Wagner , Michael B Cole , Yoon Gi Choi , Nir Yosef , Elizabeth Purdom , Sandrine Dudoit , Davide Risso , John Ngai , and Russell B Fletcher . Injury Activates Transient Olfactory Stem Cell States with Diverse Lineage Capacities . Cell stem cell , 21 ( 6 ): 775 – 790 , 12 2017 . ISSN 1875-9777 . doi: 10.1016/j.stem.2017.10.014 . URL http://www.ncbi.nlm.nih.gov/pubmed/29174333 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC5728414 . OpenUrl CrossRef PubMed Davis J McCarthy and Gordon K Smyth . Testing significance relative to a fold-change threshold is a TREAT . Bioinformatics (Oxford, England) , 25 ( 6 ): 765 – 71 , 3 2009b . ISSN 1367-4811 . doi: 10.1093/bioinformatics/btp053 . URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2654802&tool=pmcentrez&rendertype=abstract . OpenUrl CrossRef PubMed Web of Science ↵ Yongcheol Cho and Valeria Cavalli . HDAC signaling in neuronal development and axon regeneration . Current Opinion in Neurobiology , 27 : 118 – 126 , 8 2014 . ISSN 0959-4388 . doi: 10.1016/J.CONB.2014.03.008 . URL https://www.sciencedirect.com/science/article/abs/pii/S0959438814000580 . OpenUrl CrossRef PubMed ↵ Semil P. Choksi , Gilbert Lauter , Peter Swoboda , and Sudipto Roy . Switching on cilia: transcriptional networks regulating ciliogenesis . Development , 141 ( 7 ): 1427 – 1441 , 4 2014 . ISSN 1477-9129 . doi: 10.1242/dev.074666 . URL https://journals.biologists.com/dev/article/141/7/1427/46715/Switching-on-cilia-transcriptional-networks . OpenUrl Abstract / FREE Full Text Dario Magnani , Laurette Morle , Kerstin Hasenpusch-Theil , Marie Paschaki , Monique Jacoby , Stéphane Schurmans , Benedicte Durand , and Thomas Theil . The ciliogenic transcription factor Rfx3 is required for the formation of the thalamocortical tract by regulating the patterning of prethalamus and ventral telencephalon . Human Molecular Genetics , 24 ( 9 ): 2578 – 2593 , 5 2015 . ISSN 1460-2083 . doi: 10.1093/hmg/ddv021 . URL https://academic.oup.com/hmg/article-lookup/doi/10.1093/hmg/ddv021 . OpenUrl CrossRef PubMed ↵ Y Okada , S Saika , N Hashizume , S Kobata , O Yamanaka , Y Ohnishi , and E Senba . Expression of fos family and jun family proto-oncogenes during corneal epithelial wound healing . Current eye research , 15 ( 8 ): 824 – 32 , 8 1996 . ISSN 0271-3683 . doi: 10.3109/02713689609017623 . URL http://www.ncbi.nlm.nih.gov/pubmed/8921225 . OpenUrl CrossRef PubMed Web of Science Brian K. Dieckgraefe and Danielle M. Weems . Epithelial injury induces Egr-1 and Fos expression by a pathway involving protein kinase C and ERK . American Journal of Physiology-Gastrointestinal and Liver Physiology , 276 ( 2 ): G322 – G330 , 2 1999 . ISSN 0193-1857 . doi: 10.1152/ajpgi.1999.276.2.G322 . URL http://www.ncbi.nlm.nih.gov/pubmed/9950805 https://www.physiology.org/doi/10.1152/ajpgi.1999.276.2.G322 . OpenUrl CrossRef PubMed Web of Science ↵ Bradley J Goldstein , Rhea Choi , and Garrett M Goss . Multiple polycomb epigenetic regulatory proteins are active in normal and regenerating adult olfactory epithelium . Laryngoscope investigative otolaryngology , 3 ( 5 ): 337 – 344 , 10 2018 . ISSN 2378-8038 . doi: 10.1002/lio2.180 . URL http://www.ncbi.nlm.nih.gov/pubmed/30410986 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC6209616 . OpenUrl CrossRef PubMed ↵ Thomas A. Louis . Finding the Observed Information Matrix When Using the EM Algorithm . Journal of the Royal Statistical Society: Series B (Methodological) , 44 ( 2 ): 226 – 233 , 1 1982 . ISSN 00359246 . doi: 10.1111/j.2517-6161.1982.tb01203.x . URL http://doi.wiley.com/10.1111/j.2517-6161.1982.tb01203.x . OpenUrl CrossRef Web of Science ↵ Babak Alipanahi , Andrew Delong , Matthew T Weirauch , and Brendan J Frey . Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning . Nature Biotechnology , 33 ( 8 ): 831 – 838 , 8 2015 . ISSN 1087-0156 . doi: 10.1038/nbt.3300 . URL http://www.nature.com/articles/nbt.3300 . OpenUrl CrossRef PubMed Žiga Avsec , Melanie Weilert , Avanti Shrikumar , Sabrina Krueger , Amr Alexandari , Khyati Dalal , Robin Fropf , Charles McAnany , Julien Gagneur , Anshul Kundaje , and Julia Zeitlinger . Base-resolution models of transcription-factor binding reveal soft motif syntax . Nature Genetics , 53 ( 3 ): 354 – 366 , 3 2021 . ISSN 1061-4036 . doi: 10.1038/s41588-021-00782-6 . URL http://www.nature.com/articles/s41588-021-00782-6 . OpenUrl CrossRef PubMed ↵ Gemma E. Moran , Veronika Ročková , and Edward I. George . Spike-and-slab Lasso bicluster-ing . The Annals of Applied Statistics , 15 ( 1 ): 148 – 173 , 3 2021 . ISSN 1932-6157 . doi: 10.1214/20-AOAS1385 . URL https://projecteuclid.org/journals/annals-of-applied-statistics/volume-15/issue-1/Spike-and-slab-Lasso-biclustering/10.1214/20-AOAS1385.full . OpenUrl CrossRef ↵ Tianxi Li , Lihua Lei , Sharmodeep Bhattacharyya , Koen Van den Berge , Purnamrita Sarkar , Peter J. Bickel , and Elizaveta Levina . Hierarchical community detection by recursive partitioning . Journal of the American Statistical Association , pages 1 – 39 , 10 2020 . ISSN 0162-1459 . doi: 10.1080/01621459.2020.1833888 . URL https://www.tandfonline.com/doi/full/10.1080/01621459.2020.1833888 . OpenUrl CrossRef Thomas P. Minka . Estimating a Dirichlet distribution . Technical Report , 2012 . URL https://tminka.github.io/papers/dirichlet/ . View the discussion thread. Back to top Previous Next Posted March 19, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following transfactor: Transcription factor activity estimation via probabilistic gene expression deconvolution Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share transfactor: Transcription factor activity estimation via probabilistic gene expression deconvolution Koen Van den Berge , Peter Bickel , Sandrine Dudoit bioRxiv 2025.03.19.644088; doi: https://doi.org/10.1101/2025.03.19.644088 Share This Article: Copy Citation Tools transfactor: Transcription factor activity estimation via probabilistic gene expression deconvolution Koen Van den Berge , Peter Bickel , Sandrine Dudoit bioRxiv 2025.03.19.644088; doi: https://doi.org/10.1101/2025.03.19.644088 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7616) Biochemistry (17625) Bioengineering (13852) Bioinformatics (41825) Biophysics (21397) Cancer Biology (18524) Cell Biology (25417) Clinical Trials (138) Developmental Biology (13350) Ecology (19858) Epidemiology (2067) Evolutionary Biology (24277) Genetics (15581) Genomics (22459) Immunology (17698) Microbiology (40278) Molecular Biology (17134) Neuroscience (88400) Paleontology (666) Pathology (2823) Pharmacology and Toxicology (4812) Physiology (7632) Plant Biology (15106) Scientific Communication and Education (2042) Synthetic Biology (4281) Systems Biology (9807) Zoology (2266)