In recent years, a growing interest has developed around the problem of retrieving, starting from observed data, the structure of graphs representing relationships among variables of interest. In fact, reconstruction of a graphical model, known as structure learning, traces back to the beginning of the nineties, and a vast literature exists that considers the problem from various perspectives, within both frequentist and Bayesian approaches (see Drton and Maathuis (2017) for an extensive review). But a central role in the renewal of interest on structure learning equivalent to a finite mixture of singular Gaussian distributions. The authors' model, however, does not account for the count nature of the data.
From this quick tour on problems and methods, it appears evident that principled solutions to structure learning that account for the possibility of over-dispersion and/or zero-inflation are still lacking. In this paper, we try to fill this gap. We present a general framework, based on the zero-inflated negative binomial distribution, for learning the structure of a graph from single-cell RNA-seq data. We focus in particular on droplet data, for two reasons: (i) its growing popularity suggests that the majority of future studies will employ this protocol; (ii) given the high dimensionality of the problem, it is desirable to have a large sample size and droplet protocols allow the sequencing of more cells.
The remainder of this article is organized as follows. In Section 2 we introduce a motivating dataset; we describe our proposed model in Section 3 and our structure learning procedure in Section 4. One key question in the literature is whether zero inflation needs to be accounted for in the data, we offer our perspective in Section 5. After exploring the behavior of our method in simulated data in Section 6, we apply our proposal to real single-cell RNA-seq data in Section 7. Section 8 concludes the article with a discussion.
Despite the distributional challenges described in the previous section, single-cell data offer an unprecedented opportunity to discover cellular dynamics, especially in developing cell populations. Here, we study gene expression from the mouse olfactory epithelium (OE). This tissue is made of two major mature cell types, olfactory neurons and sustentacular support cells. Furthermore, a stem cell niche provides a mechanism through which the tissue is regenerated (Fletcher and others, 2017; Gadye and others, 2017) . As the aim of the study is to understand how stem cells mature into neurons following tissue damage, we focused only on the cells in the neuronal lineage. Briefly, the olfactory reserve stem cells, called Horizontal Basal Cells (HBC), become activated and subsequently develop into Globose Basal Cells (GBC) and then into immature (iOSN) and finally mature olfactory neurons (mOSN). By reconstructing the structure of the graph for each of these cell types separately, we hope to get a glimpse of the dynamic relationships between genes in neuronal development.
To this aim, we analyzed a set of 17,743 cells, assayed with 10X Genomics (v2 chemistry) after injury of the OE, to characterize HBCs and their descendants during regeneration (Brann and others, 2020) . Starting from an initial set of 25,469 cells, low-quality samples as well as potential doublets were removed as described in Brann and others (2020) . After clustering with the Leiden algorithm (Traag and others, 2019) , known marker genes were used to identify cell types. We discarded the cell types outside of the neuronal lineage (macrophages, sustentacular cells, and microvillar cells), obtaining our final dataset consisting of 7782 HBCs, 5418 activated HBCs (HBC*), 755 GBCs, 2859 iOSN, and 929 mOSN. For more details on the data preprocessing, see Brann and others (2020) .
A probabilistic graphical model requires the definition of a pair, (G, F) say. Here, G = (V, E) represents an undirected graph, where V is the set of nodes, and E = {(s, t) : s, t ∈ V, s = t} represents the set of undirected edges. Each node in the graph corresponds to a random variable X s , s ∈ V ; the existence of an edge (s, t) ∈ E indicates the dependency of the random variables X s and X t . Moreover, F represents a family of probability measures for the random vector X V , indexed by V and with support X V .
Thanks to the well known Markov properties (global, local, pairwise, see Lauritzen (1996) ), the pattern of edges in the graph translates into conditional independence properties for variables in X V , which, in turn, allow possible factorizations of F into smaller, more tractable entities. In undirected graphical models, each absent edge (s, t) in E has the role of portraying the conditional independence, X s ⊥ ⊥ X t |X V \{s,t} , and the family F is said to satisfy the pairwise Markov property with respect to G. The smallest undirected graph G with respect to which F is pairwise Markov is given the name conditional independence graph. When all variables in X V are discrete with positive joint probabilities, as is the case of this paper, the three kinds of Markov properties are equivalent, so that a factorization of the joint probability distribution with respect to the cliques (fully connected subsets of vertices) of the graph G is also guaranteed (Lauritzen, 1996, Chap. 3) .
Here, we assume that the distribution of each variable X s , conditional to all possible subsets of variables X K , K ⊆ V is a zero-inflated negative binomial (zinb) distribution: X s |x K\{s} ∼ zinb(x s ; µ s|K , θ s , π s|K |x K\{s} ), i.e., f zinb (x s ; µ s|K , θ s , π s|K |x K\{s} ) = π s|K δ 0 (x s ) + (1 − π s|K )f nb (x s , µ s|K , θ s |x K\{s} ),
(3.1)
where δ 0 (.) is the Dirac function, π s|K ∈ [0, 1] is the probability that a 0 is sampled from a distribution degenerate at zero and f nb (., µ, θ) denotes the probability mass function of the negative binomial (NB) distribution with mean µ and inverse dispersion parameter θ. We assume that ln(µ s|K ) = ν µ s|K + t∈K\{s} β µ st|K x t , (3.2) logit(1 − π s|K ) = ν π s|K + t∈K\{s} β π st|K x t . (3.3)
When both µ s|K and π s|K , are not constant, a missing edge between node s and node t corresponds to the condition β µ st|K = β µ ts|K = β π st|K = β π ts|K = 0, ∀K ⊆ V \ {s}. On the other hand, one edge between node s and node t implies that at least one of four parameters β µ st|K , β µ ts|K , β π st|K , β π ts|K is different from 0. This specification defines a family of models that includes the most common models employed with count data and embraces a variety of situations. A priori, the zinb model holds no particular rank in this family, but it is a convenient starting point from which the family of models can be established. It is quite evident that, when π s|K = 0, ∀ K ⊆ V \ {s}, the model reduces to a negative binomial distribution, which, in turn reduces to a Poisson distribution when the inverse dispersion parameter θ s tends to infinity. When π s|K > 0, zero-inflation comes into play and zero-inflated Poisson and negative binomial models can be considered. In this case, when β π st|K = 0, ∀ t ∈ K \{s}, the neighborhood of a node s is defined to be the set of effective predictors of µ s|K and consists of all nodes t for which β µ st|K = 0. On the other side, when β µ st|K = 0, ∀ t ∈ K \ {s}, the neighborhood of a node s is defined to be the set of effective predictors of π s|K and consists of all nodes t for which β π st|K = 0. In other words, the family includes models in which the structure of the graph is attributable only to one of the two parameter components, π s|K or µ s|K .
The difficulty with our model specification is that the definition of a set of conditional distributions does not guarantee the existence of a valid joint distribution, i.e., a joint distribution that possesses the specified conditionals. This might create difficulties in interpreting the resulting graph in probabilistic terms: if the joint distribution does not exist, graphical separations stored in G as a result of our model specification might not correspond to conditional independence properties on F. However, our formulation guarantees existence of the joint distribution in a number of relevant subcases. In the following theorem, we clarify conditions for existence of a joint distribution coherent with the conditional specification (see Section 1.1, Supplementary Material, for a proof).
Theorem 3.1 Let X V = (X 1 , X 2 , . . . , X p ) be a p-random vector with support X V . Assume that a set of univariate conditional probability mass functions of kind (3.1) are given for variables in X V . Then, a joint distribution having those conditionals exists if and only if θ s is constant for all s ∈ V , and all regression coefficients β µ st|K are negative, ∀K ⊆ V .
The condition on negativity of local regression coefficients in (3.2) resembles a condition known in the literature of Markov random fields known as "competitive relationship" (Besag, 1974) . Generally speaking, the presence of only negative relations among entities is quite a rare event and incapability of capturing positive dependencies might be a severe drawback in various applications. Nevertheless, the existence of a joint distribution in these specific cases assures that statistical guarantees hold for conditional approaches to structure learning such as the one used in this paper and somehow softens the hazard of the use of such algorithms outside the conditions of existence of a joint distribution.
A conditional independence graph G = (V, E) on X V can be estimated by estimating, for each node s ∈ V, its neighborhood. Hence, one can proceed by estimating the conditional distribution of X s |X V \s and fixing the neighborhood of s to be the index set of variables X N (s) on which the conditional distribution depends.
To estimate the neighborhood of each node, we employ the PC-stable algorithm, a variant of the PC algorithm first proposed by Spirtes and others (2000) . The PC algorithm starts with a complete graph on V . Marginal independencies for all pairs of nodes are tested, and edges removed when marginal independencies are found. Then, for every pair of linked nodes, independence is tested conditional to all subsets of cardinality one of the adjacency sets of the two nodes. This testing procedure is iterated, increasing in turn the size of the conditioning sets, until this reaches its maximum limit, or a limit imposed by the user. Reasons for choosing the PC algorithm are many, spanning from its consistency (assuming no latent confounders) under i.i.d. sampling (Spirtes and others, 2000) , to its ability to deal with a large number of variables and only moderately large sample sizes. The variant that we employ, PC-stable (Colombo and Maathuis, 2014) , allows to control instabilities due to the order in which the conditional independence tests are performed. To perform the tests, deviance test statistics are employed, for which a chi-squared asymptotic distribution can be obtained by standard asymptotic theory.
In what follows, let X = {x (1) , . . . , x (n) } be the collection of n samples drawn from the random vectors X V , with x (i) = (x i1 , . . . , x ip ), i = 1, . . . , n. For each K ⊆ V , let X K be the set of n samples of the |K|-random vector X K = (X t : t ∈ K), with x (i)
Starting from the complete graph, for each s and t ∈ V \{s} and for any set of variables S ⊆ {1, . . . , p}\{s, t}, we test, at some pre-specified significance level, the null hypothesis H 0 : β µ st|K = β µ ts|K = β π st|K = β π ts|K = 0, with K = S ∪ {s, t}. In other words, we test if the data support the existence of the conditional independence relation X s ⊥ ⊥ X t |X S . If the null hypothesis is not rejected, the edge (s, t) is considered to be absent from the graph. A control is operated on the cardinality of the set S of conditioning variables, which is progressively increased from 0 to p − 2 or to m, m < (p − 2).
Assume X s |x K\{s} ∼ zinb(X s ; µ s|K , θ s , π s|K |x K\{s} ), as in Equation (3.1). The conditional log-likelihood for variable X s given x K\{s} is obtained by
where ν s|K , β s|K are linked to π s|K , µ s|K through Equations (3.2) -(3.3). The estimatesν s|K ,β s|K ,θ s of the parameters ν s|K , β s|K , θ s are obtained by maximizing the conditional log-likelihood given in Equation (4.4), i.e., (ν s|K ,β s|K ,θ s ) = argmax (ν s|K ,β s|K ,θs)∈R 2|K|+1 s (ν s|K , β s|K , θ s ).
See Section 1.2, Supplementary Material, for details on the estimation procedure. A deviance test statistic for the hypothesis H 0 : β µ st|K = β π st|K = 0 can be obtained as
whereν 0 s|K ,β 0 s|K ,θ 0 s are the maximum likelihood estimates of the parameters under H 0 . It is readily available that D s|K is asymptotically chi-squared distributed with 2-degrees of freedom under the null hypothesis, provided that some general regularity conditions hold.
Remark 1. On assuming faithfulness of the node conditional distributions to the graph G, consistency of the algorithm can be proved in the case of competitive relationships in µ s |K by suitably modifying results in . We recall that a distribution P X is said to be faithful to the graph G if for all disjoint vertex sets A, B, C ⊂ V it holds
where A ⊥ ⊥ G B|C means that A and B are separated in G by C. Thanks to the equivalence between local and global Markov properties, faithfulness of the local distributions guarantees faithfulness of the joint distribution.
Remark 2. A large sample size, as typical in the applications at hand, impacts on the actual significance level of individual tests. Moreover, a multiplicity of tests are performed by the algorithm. For this reason, we advice to set the nominal level of the test α to α n = 2(1 − Φ(n b )), where 0 < b < 1/2 is related to the average neighborhood size. This choice is based on results in and guarantees that the probability that a type I or II error occurs in the whole testing procedure goes to zero as n → ∞, i.e., it asymptotically controls the family-wise error rate of all potential tests that could be done.
Remark 3. The chosen learning strategy has some advantages over alternative approaches based on sparse regressions (see also for an extended discussion in the Poisson case). Sparsity can be easily implemented by a control on the conditional set size, instead of a control on parameter magnitudes, which can lead to over-shrinkage. Moreover, it offers computational advantages, especially when sparse networks are the target of inference. Finally, hypothesis testing is scale-invariant, i.e., is not affected by scale transformations of regressors.
The need for modeling zero-inflation in single cell data is a question at the core of an ongoing debate, with several authors arguing that the negative binomial distribution is sufficient to fit single-cell RNA-seq data when unique molecular identifiers are used (Vieth and others, 2017; Townes and others, 2019; Svensson, 2020). Indeed, the ability to distinguish between a non zero-inflated distribution and zero-inflated alternatives highly depends on the relative size of the parameters of the distributions.
To gain a better understanding of this problem, we have tried to assess the misspecification cost due to assuming a non zero-inflated distribution when zero inflation is in action. To this aim, we confined ourselves to a univariate case with no covariates, fixed a non zero-inflated model and measured the model misspecification cost occurring when using its zero-inflated counterpart by using the squared Hellinger distance as loss function. Such a loss function should, in principle, indicate, in an inferential sense, how far apart the two distributions are.
To this aim, let Y = {0, 1, 2, . . . , +∞} be the support of a discrete variable Y . We consider for Y a true probability distribution P (y; φ 0 ), φ 0 ∈ Φ, as well as a family F = {Q(y; ψ), ψ ∈ Ψ}, Φ ⊆ Ψ, of zero-inflated versions of the true probability distribution P (y; φ 0 ). The squared Hellinger distance between two probability distributions P and Q is defined as
where p y = P (y; φ 0 ) and q y = Q(y; ψ).
In particular, assume that P is a NB distribution with φ 0 = (µ 0 , θ 0 ) and Q is a zinb distribution, defined as Q(y; ψ) = πδ 0 (0) + (1 − π)P (y; φ 0 ), with ψ = (φ 0 , π). Hence, the squared Hellinger distance of P and Q can be written as Table 1 shows the value of the Hellinger distance in a number of cases. As expected, the distance increases with the probability of zero inflation π. However, when the inverse dispersion parameter θ 0 and/or the mean µ 0 are small, the distance between the distributions is small even in the case of moderate to large π. In fact, when µ 0 and θ 0 are both small (low mean and high variance) the two distributions are close even when π = 0.9.
As, broadly speaking, maximum likelihood estimators and minimum Hellinger distance estimators are asymptotically equivalent, it emerges that, in inferential terms, the degree of zero inflation of a true model could be difficult to ascertain, as suitable choices of the parameters of the non contaminated component may possibly absorb the excess of zeros generated by the contamination. This is despite identifiability of the zinb model (see Section 1.3, Supplementary Material, for a proof). These remarks might contribute to the ongoing debate about existence of zero-inflation from a novel perspective.
We devote this section to the empirical study of consistency of the proposed algorithms. In particular, we concentrate on the ability of proposed methods to recover the true structure of the graphs. We also list the running time of each algorithm. As measures of the test's accuracy, we adopt three criteria including Precision P ; Recall R; and their harmonic mean, known as F 1 -score, respectively defined as
where TP (true positive), FP (false positive), and FN (false negative) refer to the number of inferred edges (Liu and others, 2010) . The considered algorithms are listed below, along with specifications, if needed, of tuning parameters. For all PC-like algorithms, we let the maximum cardinality of conditional independence set be m = 8 for p = 10 and m = 3 for p = 100.
-PC-zinb1: zinb models in which the structure of the graph is attributable to both of the two parameter components µ s|K and π s|K ;
-PC-zinb0: zinb models in which the structure of the graph is attributable to only the parameter component µ s|K and consider π s|K as a constant (i.e., β π st|K = 0, ∀ t ∈ K \ {s}, ∀ s ∈ V ); -PC-nb: Negative binomial model, i.e., the special case of zinb models where π s|K = 0; -PC-pois: Poisson model .
For two different cardinalities (p = 10 and p = 100), we consider three graphs of different structure: (i) a scale-free graph, in which the node degree distribution follows a power law; (ii) a hub graph, where each node is connected to one of the hub nodes; (iii) a random graph, where the presence of the edges is drawn from independent and identically distributed Bernoulli random variables.
To construct the scale-free and random networks, we employed the R package igraph (Csardi and others, 2006) . For the scale-free networks, we followed the Barabasi-Albert model with constant out-degree of the vertices ν = 2 for p = 10 and ν = 0.2 for p = 100. For the random networks, we followed the Erdos-Renyi model with probability to draw one edge between two vertices γ = 0.3 for p = 10 and γ = 0.03 for p = 100. To construct the hub networks, we assumed 2 hub nodes for p = 10, and 5 hub nodes for p = 100. See Supplementary Figure S1 and Supplementary Figure S2 for representative plots of the three chosen graphs for p = 10 and p = 100, respectively.
For the given graphs, 50 datasets were sampled with four different sample sizes, n = {100, 200, 500, 1000} for p = 10, and three different sample sizes, n = {200, 500, 1000} for p = 100. To generate the data, we followed the approach of the Poisson models in . Let X ∈ R n×p be the set of n independent observations of random vector X. Then, X is obtained from the following model X = YA + , where Y = (y st ) is an n × (p + p(p − 1)/2) matrix whose entries y st are realizations of independent random variables Y st ∼ zinb(µ, θ, π) (or nbinom(µ, θ); or Pois(µ)) and = (e st ) is an n × p matrix with entries e st which are realizations of random variables E st ∼ zinb(µ nois , θ, π) (or nbinom(µ nois , θ); or Pois(µ nois )). This approach leverages the additive property of these distributions and allows us to generate the required dependencies. In particular, let B be the adjacency matrix of a given true graph, then A takes the following form A = [I p ; P (1 p tri(B) T )] T . Here, P is a p × (p(p − 1)/2) pairwise permutation matrix, denotes the elementwise product, and tri(B) is the (p(p − 1)/2) × 1 vectorized upper triangular part of B (Allen and Liu, 2013).
Figures 1 and 2 show the Monte Carlo means of the F 1 -scores for each of the considered methods with p = 100 and low signal-to-noise ratio (µ noise = 0.5), at high (µ = 5) and low (µ = 0.5) mean levels, respectively. Each value is computed as the average of the 50 values obtained by simulating 50 samples for the model corresponding to each network. Monte Carlo means of Precision P , Recall R, and F 1 -score are given in Tables S1-S3.
The two values of µ = {5, 0.5} were chosen to mimic typical values observed in real full-length and droplet-based datasets, respectively. In fact, the mean expression level of transcription factors in the dataset presented in Section 2 is 0.67 (median 0.14), while the mean expression level of transcription factors in a similar experiment performed with a full-length protocol (Fletcher and others, 2017) is 32.03 (median 7.89).
These results indicate that the PC-zinb1 algorithm and its variants (PC-zinb0, PC-nb, PCpois) are consistent in terms of reconstructing the structure from given data. In fact, when the model is correctly specified, the F 1 -scores of the algorithms are close to 1 when n 1000 in all scenarios. This means that the proposed algorithm is able to recover the underlying graph from the given data for both low (Fig. 2 ) and high ( Fig. 1 ) mean levels.
When the data are generated with a high mean level (µ = 5), the PC-pois algorithm performs well only when it is the true model, i.e., for data generated from Poisson random variables (Fig.  1; Supplementary Table S3 ). In the other scenarios, PC-pois often shows a low Precision ( Fig.  1 ; Supplementary Tables S1 and S2). This result is expected since the node conditional Poisson distributions are unable to model the over-dispersion generated by the (zero-inflated) negative binomial distributions.
On the other end of the spectrum, the more general zinb models work well in all scenarios ( Fig. 1 ; Supplementary Tables S1 -S3). This is not surprising as the data are generated according to models (e.g, Poisson, NB) that can be seen as special cases of the zinb distribution, which means that in all tested scenarios the zinb model is correctly specified.
The PC-nb algorithm, based on the negative binomial assumption, performs reasonably well ( Fig. 1 ; Supplementary Tables S1 -S3). However, in the hub graph (middle row of Fig. 1 ), its performances are slightly worse than the zinb models, showing low Precision when the true data generating distribution is node conditional zinb ( Fig. 1; S1 ). This result indicates that a zero inflated negative binomial model may be needed when the mean is large (Risso and others, 2018) .
As we expected from the considerations reported in Section 5, the performances of the variants of PC-zinb are quite similar to each other when the mean and the dispersion parameter are both small, i.e., when the data are characterized by low mean and high variance (µ = 0.5, θ = 0.5; Fig.  2 ; Supplementary Tables S1 -S3). This might be explained by the fact that a suitable choice of the parameters may allow non-zero inflated models to absorb the excess of zeros (see Section 5 for more details). Therefore, when applying our approach on this type of data, one should use the simplest variant, (i.e., PC-pois) to leverage the better computational performance (see last column of S1 -S3).
Finally, we see no difference in the performance of the PC-zinb variants (PC-zinb1 and PC-zinb0). This is perhaps not surprising, as we simulated the same structure of the graph for both µ and π. These results suggest that the information inferred from µ is sufficient to reconstruct the correct graph in this case.
We have focused here on p = 100, as this setting is closer to our real application. The results for p = 10 are reported in Supplementary Figures S1 and S2 and Supplementary Tables S4-S6 and lead to similar conclusions.
We demonstrate our method on the motivating example dataset described in Section 2. Our analysis focuses on the total set of 1543 known transcription factors in mouse, which are thought to regulate the observed differentiation processes. We furthermore focus on the differentiation path starting at the activated HBC (HBC*) stage (i.e., activated stem cells upon injury) up to mature neurons, therefore investigating the entire neuronal lineage in the trajectory of this dataset. As previously discussed in Section 2, we expect four different cell types along this path, being respectively HBC*, GBC, iOSN and mOSN, and we estimate the structure of the graph for each of these cell types.
The average degree of the graphs is highest at the activated stem cell stage, with an average degree of 4, and decreases as cells develop to mature neurons, with average degrees of 3.9, 3.3 and 3.5 for the GBC, iOSN and mOSN networks, respectively. To interpret the graph structure, we focus on the 2-core of each network, i.e., we retain TFs that are associated with at least two other TFs, a preprocessing step that helps in understanding the core structure (Wang and Rohe, 2016).
We identify communities in each graph using the Leiden algorithm (Traag and others, 2019) and, in order to validate the associations discovered by PC-zinb, we interpret each of the communities by computing overlap with known functional gene sets in the MSigDB database (Subramanian and others, 2005; Liberzon and others, 2015), see Supplementary Material Section 2 for details. The interpretation of these communities relies on known processes involved in the development of the olfactory epithelium as found by previous research (e.g., Fletcher and others (2017); Gadye and others (2017)).
In the HBC* cell type cluster, cells have been injured ∼ 24h ago, so we expect response to injury, and stem cells actively preparing for differentiation, as well as replication to produce more stem cells to repair the epithelium. Four communities are discovered in the association network ( Figure 3 ), broadly involved in either cell cycle, epigenetic mechanisms and (epithelial) cell differentiation (Supplementary Table S8 ). These communities reflect the need to divide in order to produce more cells, epigenetic mechanisms that are likely required to activate molecular processes upon injury, and the differentiation of stem cells to restore the damaged epithelium.
In the GBC cell type, we expect cells to proliferate to produce immature neurons. We discover four communities (Figure 3 ), broadly involved in DNA replication, cell proliferation, signaling, expression regulation and cell differentiation (Supplementary Table S9 ). Relevant pathways, such as the P53 and notch signaling pathways, are also recovered for specific communities, and have previously been found to be involved in neurogenesis in neuroepithelial stem cells (Marin Navarro and others, 2020; Wang and others, 2011).
In the immature olfactory sensory neuron (iOSN) stage, we expect basal cells to start developing into immature neurons. Four communities are discovered (Figure 3 ), of which one community comprises the majority of the graph, i.e., 63% of all TFs retained in the graph, and importantly is involved in neurogenesis (Supplementary Table S10 ). Other, smaller, communities are enriched in processes such as cell and axon growth, wound healing, signaling and cell population maintenance.
Finally, in the mature olfactory sensory neuron stage (mOSN), we expect the final differentiation to functional neurons. Five communities are discovered (Figure 3) , again with very different sizes. The largest communities are enriched in broader processes related to chromatin organization and transcription, possibly reflecting the basic changes required for cells to develop into and maintain at the mature stage (Supplementary Table S11 ). The third largest community is enriched specifically in the TGF-Beta pathway, known to be required for neurogenesis, and to modulate inflammatory responses (Meyers and Kessler, 2017) .
Taken together, these results confirm previously known processes associated with differentiation of HBCs into mature neurons upon injury, with relevant processes highlighted by communities of transcription factors. Furthermore, while the community detection results are useful to validate the estimated graphs, they also provide a gateway to more detailed analysis, e.g., investigation of hub genes (e.g., Chen and others (2018)) or master regulators of development (e.g., Sikdar and Datta (2017) ), therefore unlocking powerful interpretation of single-cell RNA-seq datasets.
In this work, we have introduced PC-zinb, a class of constraint-based algorithms for structure learning, supporting possibly overdispersed and zero-inflated count data. In focusing on these two nonstandard but realistic situations, our framework goes beyond what has so far been proposed in the literature. Moreover, by leveraging the proposal in Nguyen and Chiogna (2018) -shown to be competitive with state-of-the-art methods supporting count data -we inherit the benefits of that approach, most notably: the existence of a theoretical proof of convergence of the algorithm under suitable assumptions; an easy implementation of sparsity by a control on the number of variables in the conditional sets; invariance to feature scaling. On the synthetic datasets considered in Section 6, we showed that the algorithms work well in terms of reconstructing the structure from given data for large enough sample sizes, while providing biologically coherent information and insight on the real dataset analyzed in Section 7.
Our simulation studies allow us to derive various recommendations on the use of PC-zinb. Clearly, these do not rule out sensitivity analyses with respect to both model specification and tuning of the algorithms, which remain an important part of the model criticism process. A control of the level of significance of the tests with respect to the sample size, n, and the expected size of the neighborhood of each node, b, is highly recommended to guarantee good reconstruction abilities. As in real applications knowledge of the expected size of the neighborhood might be difficult to elicit, it may be prudent to try a range of values for b, and check stability of results. This might also generate a sequence of models of decreasing complexity for increasing values of b and whose dynamic might also point researchers to the most significant connections.
If only the structure of the graph is of interest, irrespective of the strength of the links among variables, we suggest making use of the Poisson variant of the algorithm when the mean of the variables is small, so as to reduce computational complexity. Moreover, when the mean of the variables is small, presence of zero-inflation might not influence reconstruction abilities of the algorithms, as also confirmed by the small study on zero-inflation in Section 5. In these situations, we recommend using, at least in the first instance, non zero-inflated models.
If the null hypothesis H 0 : β µ st|K = β π st|K = 0 fails to be rejected, PC-zinb will remove the edge between variables s and t. While such a procedure can only be justified in settings with high power, our simulation study shows that, even in settings with small sample sizes, our algorithm is able to achieve high power, and the correct underlying structure of the graph can be learned successfully.
Our real data analysis, aimed at assessing biological validity of the reconstructed network, has demonstrated the great importance of finding meaningful visualizations of large complex networks. Our proposal, based on a search for communities of variables and their association to gene ontologies via enrichment analysis, allowed us to confirm both biological interpretability of the estimated structure, and to contribute to our understanding of which and where biological processes are occurring.
The methods presented in this article are available in the learn2count R package, available at https://github.com/drisso/learn2count. The code to reproduce the analyses of this paper is available at https://github.com/drisso/structure_learning.
Supplementary material includes proofs and mathematical details, details on the real data analysis, and supplementary figures and tables. Yang, E., Ravikumar, P. K., Allen, G. I. and . On poisson graphical models.
In: Advances in neural information processing systems. pp. 1718-1726.
Yang, E., Ravikumar, P. K., Allen, G. I. and Liu, Z. 
In Section 1.1, we provide the proof of Theorem 3.1 (Section 3 of main paper); Section 1.2 describes the maximum likelihood estimation procedure used in zinb regressions; in Section 1.3 we prove the identifiability of zinb models (Section 3 of main paper).
In what follows, we will first consider the case K = V . Results for general case K ⊆ V will be then deduced. We note that the set of explanatory variables X K\{s} in the generalized linear model X s given X K\{s} does not include variable X t , with t ∈ V \K. Suppose we zero-pad the true parameter β s|K ∈ R 2|K|−2 to include zero weights over V \K, then the resulting parameter would lie in R 2|p−1| . To simplify notation, we write (µ j , π j , β µ jk , β π jk ) instead of (µ j|V , π j|V , β µ jk|V , β π jk|V ).
Proof. It is easy to show that for x = (x 1 , x 2 , . . . , x p ) and y = (y 1 , y 2 , . . . , y p ) one can write P (x) P (y) = p j=1 P (x j |x 1 , . . . , x j−1 , y j+1 , . . . , y p ) P (y j |x 1 , . . . , x j−1 , y j+1 , . . . , y p ) .
1 arXiv:2011.12044v1 [q-bio.GN] 24 Nov 2020
Let y = (0, 0, . . . , 0). It holds
Consider the following two cases.
.
To prove existence of the joint distribution, it is sufficient to control the term
Indeed, taking the sum over x 1 , x 2 , . . . , x p , one gets x1,...,xp
for some q < 1. Therefore, when x j is large enough, from Equation 1.2, one has x1,...,xp
that proves existence of the joint distribution.
• At least one parameter β µ jk > 0.
Without loss of generality, one can assume β µ p(p−1) > 0. Then, µ p−1 p → ∞, when x p−1 → ∞. From Equation 1.1, one has
Taking the sum over x 1 , x 2 , . . . , x p of the term
(1.3)
-If β π p(p−1) > 0, then 1 − π p−1 p = exp{β π p + β π p(p−1) x p−1 } 1 + exp{β π p + β π p(p−1) x p−1 } → 1 when x p−1 → ∞. 
Therefore, the joint distribution exists if and only if all regression coefficients β µ st|K are negative, ∀K ⊆ V .
Given a set of data X, for each s ∈ K ⊆ V , we estimate parameters ν s|K , β s|K , θ s by maximizing the local log-likelihood function
To solve this standard regression problem, we alternate the estimation of the dispersion parameter θ s and of parameters (ν s|K , β s|K ), as done in Risso and others (2018) . Specifically, after initialization, we iterate the following two steps until convergence (for the reader's convenience, we drop the index K). Note that this procedure can be performed independently and in parallel for each regression performed within the structure learning strategy.
Step 1: dispersion optimization. Let ln θ s = ζ s , we estimate parameter ζ s by solving the following problemζ
A first initial value for ζ s , ζ 0 s say, is found by applying the method of moments. If the moment estimator does not exist, as in the Poisson case, ζ 0 s is assigned to the value that maximizes the objective function s (ν s ,β s , ζ s ) over a suitable interval. In practice, we use optimize function in R over the range [-100, 100] . Once the initial value ζ 0 s is fixed, we optimize s (ν s ,β s , ζ s ) by a quasi-Newton optimization scheme.
The gradient of the objective function used by the optimization procedure is found as follows. The derivative of the negative binomial log density is:
where Ψ(θ) = Γ (θ)/(Γ(θ)) is the digamma function. Hence, the derivative of the zinb model for π ∈ [0, 1] can be written as follows:
• for y > 0, f zinb (y, µ, θ, π) = (1 − π)f nb (y, µ, θ), so that
(1−π)θ θ + 1 , Therefore, the derivative of the objective function w.r.t. ζ s for each response X s , s = 1, . . . , p, is
Step 2: log-likelihood optimization. Given a collection of n samples drawn from the random vector X, X say, and dispersion parameterζ s ∈ R, we estimate parameter (β s , ν s ) by solving the following problem
For convenience, we reparameterize the log-likelihood as follows, as done in Risso and others (2018) :
Therefore, Problem (1.4) translates into finding a set of vectors (a µ , a π ) that locally maximizes the log-likelihood of a zinb model for a vector of count X s , i.e., F (a µ , a π ) = s (ν s , β s ,ζ s ).
Starting from an initial value, we perform a local maximization of this function F using BFGS method (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970) . We first explicit the derivatives of the log-likelihood w.r.t. parameters (µ, π) in the zinb distribution.
We note that f zinb (y; µ, θ, π) = πδ 0 (y) + (1 − π)f nb (y, µ, θ),
Depending on whether y is equal to zero, the above given expressions become:
-y > 0: here, δ 0 (y) = 0 and f zinb (y, µ, θ, π) = (1 − π)f nb (y, µ, θ), so that one has ∂ ∂µ ln f zinb (y, µ, θ, π) = y θ − y + θ θ + µ ∂ ∂π ln f zinb (y, µ, θ, π) = −1 1 − π .
-y = 0: here, δ 0 (y) = 1 and f zinb (0, µ, θ, π) = π + (1 − π)f nb (0, µ, θ), and one has
Using standard calculus for the differential of compositions and the fact that:
one obtains
where G and H are the n-dimensional vectors given by
Remark 1.1 Note that for each s ∈ V , the dispersion parameter θ s is constant, i.e., θ s does not depend on conditional sets K. Hence, in practice, we estimate θ s by applying the above given procedure for K = V . Then, we use this estimateθ s as the input of the second step, log-likelihood optimization.
Remark 1.2 It is worth to note that: (i) for negative binomial models, i.e., π s|K = 0, ∀s ∈ K ⊆ V , parameters are estimated similarly by considering a π as a vector 0, and A π as a matrix 0; (ii) negative binomial distributions belong to the exponential family. Hence, an alternative way to estimate parameters is the iteratively reweighted least squares (IRWLS) algorithm. This approach does not depend on the choice of initial values for parameters. As a result, PC-nb using IRWLS algorithm, called PC-nbglm could be more accurate in some scenarios with significantly reduced computational time (see Table S7 ). In practice, when dealing with high dimensional data, we suggest users to use the PC-nbglm algorithm.
In this section, we prove that the zinb model is identifiable.
Theorem 1.3 Assume X follows a zinb distribution, i.e., X ∼ zinb(π, θ, µ). Then, the density function of X is identifiable.
Proof. We prove the theorem by contradiction. Assume there exists (π 1 , θ 1 , µ 1 ) = (π 2 , θ 2 , µ 2 ) such that X 1 ∼ zinb(π 1 , θ 1 , µ 1 ) with density function f 1 (., π 1 , θ 1 , µ 1 ); X 2 ∼ zinb(π 2 , θ 2 , µ 2 ) with density function f 2 (., π 2 , θ 2 , µ 2 ), and f 1 (., π 1 , θ 1 , µ 1 ) ≡ f 2 (., π 2 , θ 2 , µ 2 ). For s = 1, 2, one has f s (x s , µ s , θ s , π s ) = π s δ 0 (x s ) + (1 − π s )f nb (x s , µ s , θ s ). Let x s = 0, 1, 2, 3, then
One has f 1 (2, π 1 , θ 1 , µ 1 ) f 1 (1, π 1 , θ 1 , µ 1 ) = f 2 (2, π 2 , θ 2 , µ 2 ) f 2 (1, π 2 , θ 2 , µ 2 ) ⇒ Γ(2 + θ 1 ) Γ(1 + θ 1 )
(1.6) As Γ(z + 1) = zΓ(z), Equation 1.6 yields
(1.7)
Similarly,
(1.8)
From Equation (1.7) and Equation (1.8), one gets
(1.9)
Finally, Equation (1.5), (1.7), and (1.9) yield (π 1 , θ 1 , µ 1 ) = (π 2 , θ 2 , µ 2 ), that contradicts the Assumption.
Our analysis focusses on the total set of 1543 transcription factors (TF), which are known to play a central role in the regulation of other genes (Fletcher and others, 2017) . As measurements were zero-inflated and highly skewed, with total count volumes depending on experimental condition, standard preprocessing was applied to the data as in ; . In particular, we filtered out the genes with the mean smaller than 0.005; normalized the data by 95% quantile matching; selected the top 1000 cells with the highest mean from the groups that had more than 1000 cells; and used a power transform X α (with α = 0.5). Normalized data were used as input to PC-zinb1 with a significance level of α = 2(1−Φ(n 0.15 )) , and the maximum cardinality of conditional independence set m = 3. Once estimated, we extracted the 2-core of the networks by iteratively excluding transcription factors with fewer than 2 edges, leaving a network with 926, 991, 465 and 772 genes for the HBC*, GBC, iOSN and mOSN cell types, respectively. For each network, we performed community detection using the Leiden algorithm (Traag and others, 2019), where we supervise the resolution parameter to be 0.4 for HBC*, GBC and mOSN, and 0.25 for the iOSN network, as a resolution parameter of 0.4 resulted in a too high number of communities there. The transcription factors in each of these communities were then used in an overrepresentation analysis using the biological processes hallmark gene sets of the MSigDB database Liberzon and others, 2015) , where only gene sets containing at least 5 genes were retained. We ranked gene sets using a Fisher's exact test, and only the top 10 enriched gene sets were interpreted. Fig. S3 . F1-score of the considered algorithms for the three types of graphs in Supplementary Figure S1 with p = 10, µ = 5, θ = 0.5: random (top); hub (middle); scale-free (bottom). PC-zinb1: zinb model in which the structure of the graph is attributable to both of the two parameter components µ s|K and π s|K ; PC-zinb0: zinb model in which the structure of the graph is attributable to only the parameter component µ s|K and consider π s|K as a constant; PC-nb: Negative binomial model, i.e., the special case of zinb models where π s|K = 0; PC-pois: Poisson model of . Figure S2 (p = 100), under NB node conditional distribution mean µ = 0.5, 5; θ = 0.5 and levels of noise µnoise = 0.5. The levels of significance of tests α = 2(1 − Φ(n 0.2 )) for n = 500, 1000, and α = 2(1 − Φ(n 0.225 )) for n = 200. Figure S2 (p = 100), under Poisson node conditional distribution with mean λ = 0.5, 5; and levels of noise λnoise = 0.5. The levels of significance of tests The levels of significance of tests α = 2(1 − Φ(n 0.2 )) for n = 500, 1000, and α = 2(1 − Φ(n 0.225 )).
PC-zinb0 Figure S1 (p = 10), under zinb node conditional models with mean µ = 0.5, 5; θ = 0.5, and levels of noise µnoise = 0.5. The levels of significance of tests α = 2(1 − Φ(n 0.15 )). Figure S1 (p = 10), under NB node conditional models with mean µ = 0.5, 5; θ = 0.5, and levels of noise µnoise = 0.5. The levels of significance of tests α = 2(1 − Φ(n 0.15 )). Figure S1 (p = 10), under Poisson node conditional model with mean λ = 0.5, 5; and levels of noise λnoise = 0.5.
The levels of significance of tests α = 2(1 − Φ(n 0.15 )). Supplementary Table S7 . Monte Carlo means of Precision (P), Recall (R), F1-score, and runtime obtained by simulating 50 samples from three different types of graphs in Figure S2 (p = 100), under zinb and NB node conditional distributions mean µ = 5; θ = 0.5; π = 0.7 and levels of noise µnoise = 0.5. The levels of significance of tests The levels of significance of tests α = 2(1 − Φ(n 0.2 )) for n = 500, 1000, and α = 2(1 − Φ(n 0.225 )) for n = 200. 