Estimating graphical models for count data with applications to single-cell gene network
Authors: Feiyi Xiao
Junjie Tang
Huaying Fang
Ruibin Xi
Venue: NeurIPS 2022 Accept (2022)
Keywords: Convergence rate
Graphical model
Network inference
Poisson log-normal model
Probabilistic Methods
Single-cell RNA-Seq.

Abstract: Graphical models such as Gaussian graphical models have been widely applied for direct interaction inference in many different areas. In many modern applications, such as single-cell RNA sequencing (scRNA-seq) studies, the observed data are counts and often contain many small counts.  Traditional graphical models for continuous data are inappropriate for network inference of count data. We consider the Poisson log-normal (PLN) graphical model for count data and the precision matrix of the latent normal distribution represents the network. We propose a two-step method PLNet to estimate the precision matrix. PLNet first estimates the latent covariance matrix using the maximum marginal likelihood estimator (MMLE) and then estimates the precision matrix by minimizing the lasso-penalized D-trace loss function. We establish the convergence rate of the MMLE of the covariance matrix and further establish the convergence rate and the sign consistency of the proposed PLNet estimator of the precision matrix in the high dimensional setting. Importantly, although the PLN model is not sub-Gaussian, we show that the PLNet estimator is consistent even if the model dimension goes to infinity exponentially as the sample size increases. The performance of PLNet is evaluated and compared with available methods using simulation and gene regulatory network analysis of real scRNA-seq data.

Estimating graphical models for count data with applications to single-cell gene network
Feiyi Xiao
Junjie Tang"
School of Mathematical Sciences
School of Mathematical Sciences
Peking University xiaofeiy11217$pku edu. cn
Peking University
junjie tangâpku edu.
Huaying Fang
Beijing Advanced Innovation Center for Imaging Theory and Technology Capital Normal University
huayingfang@ho com
Ruibin
School of Mathematical Sciences
Center for Statistical Science Peking University
ruibinxi@math pku. edu. cn
Abstract
Graphical models such as Gaussian graphical models have been widely applied for direct interaction inference in many different areas. In many modern applications, such as single-cell RNA sequencing (scRNA-seq) studies, the observed data are counts and often contain many small counts. Traditional graphical models for continuous data are inappropriate for network inference of count data. We consider the Poisson log-non log-normal (PLN) graphical model for count data and the precision matrix of the latent normal distribution represents the network. We propose a two- step method PLNet to estimate the precision matrix. PLNet first estimates the latent covariance matrix using the maximum marginal likelihood estimator (MMLE) and then estimates the precision matrix by minimizing the lasso-penalized D-trace loss function. We establish the convergence rate of the MMLE of the covariance matrix and further establish the convergence rate and the sign consistency of the proposed PLNet estimator of the precision matrix in the high dimensional setting. Importantly, although the PLN model is not Gab--aassian, we show that the PLNet estimator is consistent even if the model dimension goes to infinity exponentially as the sample size increases. The performance PLNet evaluated and compared with available methods using simulation and gene regulatory network analysis of real scRNA-seq data.
Introduction
Gaussian graphical model (GGM) 15. 28
is widely used for understanding the complex interactions
of the observed variables various fields [26] 14 (33]. GGM assumes that each sample is drawn from multivariate Gaussian distribution, in which the precision matrix (i.c., the inverse of the covariance matrix) represents the network. In GGM, two nodes are connected if the corresponding
'Equal contribution. Corresponding author: Ruibin Xi (ruibinxi math. mk....pu..ddu.c)
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

element in the precision matrix is non-zero. Based on the sparsity assumption of the precision matrix, researchers have proposed many methods for estimating the precision matrix of the GGM and established consistency theories, such as methods by maximizing the penalized log-likelihood [28 by solving an equivalent regression problem with the lasso penalty [5] 18 or by minimizing
smooth convex loss function (called D-trace loss) with the lasso penalty (3) Due to its high interpretability, well-established theoretical properties, and computational advantages, GGM has been adopted in many applications.
This paper is motivated by the gene regulatory network analysis based on the single-cell RNA sequencing (scRNA-seg) data. The recent breakthrough of scRNA-seq technologies has provided tremendous opportunities for understanding transcriptional states and activities and for uncovering gene regulatory networks at single cell level. However, gene expressions from scRNA-seq are count data and the Gaussian assumption is inappropriate. In particular, for the recent unique molecular identifier (UMI) based scRNA-seq data [ 32]. the expression counts contain many zeros. Trans- formations such as taking logarithms cannot make the Gaussian model a good approximation to the distribution of scRNA-seq expression data and may distort the correlation structure when there are many zeros in scRNA-seq data.
Poisson distribution is a natural choice for modeling count data. Researchers have generalized the univariate Poisson distribution to multivariate distributions for network analysis [9]. Some such models are the Poisson graphical model [2] and its variants generalized Poisson graphical models [25]. But these generalized Poisson graphical models can not handle the over-dispersion often observed in real data. Negative binomial distributions are often used to account for the over-dispersion [19], but it is more difficult to generalize negative binomial distributions to describe the network structure in the multivariate count data. Here, we use the multivariate Poisson normnolaal (PLN) distribution for graphical modeling of count data. The PLN distribution is a hicrarchical model of Poisson and multivariate log-normal distributions Q. A random vector Y (Yi,..., ,,,...,,)''' is from PLN distribution, if conditional on latent vector X = (X1, ... X, Xp)² with (log(Xi X1)..... log( X,,) x N(H, E), elements of Y are independent Poisson random variables with mean parameters XI,..., Xp-
The major advantage of the PLN model is that, similar to the GGM, the precision matrix é E-1 can represent the network. Network recovery can be achieved by estimating O. In addition. PLN distributions allow over-dispersion and hence are suitable for modeling count data with over- dispersion. In gene regulatory network analysis of scRNA- -sqqA--eeqaa data, the PLN graphical model has clear biological explanation. Previous experimental researches showed that gene expressions of single cells follow normnol distributions 4 The latent variables XI...., Xp. representing the true expressions of p genes in a single cell, can be reasonably assumed to be jointly oorrolllly distributed. The precision matrix of the log-normal distribution represents gene-gene interactions in single cells. The observed variables Y are measurements of the underlying true expressions X. Largely speaking, the log-normal layer of the PLNN model captures the biological fluctuation of gene expressions and the Poisson layer accounts for the technical and measurement noises. Only the biological fluctuation reflects gene-gene interactions and the regulatory network is the precision matrix of the latent log-normal model.
few algorithms have been developed for estimating the precision matrix of the PLN model in high- dimensional settings. Compared with GGM, the likelihood of the PLN model is more complicated since it involves a high- nigh-dimensional integration and does not have closed form. Maximizing the penalized log-likelihood of the PLN model is very difficult. Wu et al. 24 used Monte Carlo methods to approximate the log- likelihoodd and estimate the precision matrix by maximizing the penalized approximated log- likeliiekoood of the PLNN model. Chiquet et al. [5] developed a computationally more appealing method based on the variational approximation. However, these methods are all based on approximations of the logeliiel and the accuracies of these approximations need to be further elaborated. More importantly, no convergence theory has been developed for these precision matrix estimators.
In this article, we propose a two-step approach named PLNet that first estimates the covariance matrix of the latent log- normal variables in the PLN model using the maximum marginal likeli- hood estimator (MMLE) and then estimates the sparse precision matrix by minimizing the lasso penalized D-trace loss 30 One advantage of this two-step approach is that it avoids computing high-dimensional numerical integration that is involved in the log-likelihood of the multivariate PLNN model. Minimizing the penalized D-trace loss is computationally cost-effective. Thus, this estimator
2

---

is generally computationally more efficient. Furthermore, we show that, under mild conditions, the estimator given by PLNet consistent estimator of the precision matrix high dimensional settings. Comprehensive simulation studies show that PI.Net provides more accurate estimates and computationally more efficient than available methods. We also demonstrate the application of PLNet to seRNA-seq data.
2 Methods
2.1 Notation
For column vector a with the i-th entry as ||al| = (a"a)²/2 and |a|l; = = lasl. For matrix A with the (i.j)-th entry A.j. ||A|| = max,,j |Ay] be the element-wise 4.0 -norm. ||A||| Σng |4i
1/2 be the Frobenius norm, ||A||| || All// be the operator norm, ||-||] Σies/ |Aspl. || =
be the nnnrrmm || = max; Eg |Ay)) be the norm,, ||A||= = [iis |Ais)
Enj #A) be the number of nonzero entries where 1 (-) is the indicator function, tr (A) = E; Asi be the trace, and A = 0 means A is positive semi-definite.
2.2 The PLN graphical model
Let Y = [Kyll be the observed count matrix, Y, = -Y Yip)" be the observed count data of the ith sample, and X, (Xil,... Xip) be the latent random vector. In scRNA-seq data, Y, and X,j are the observed expression and the underlying "true" expression of the jth gene in the ith cell, respectively. We assume that conditional on X;, elements of Y; are independent Poisson random variables with mean S,X,, (1 S S p). where S; is a scaling factor sampled from one distribution g(S) and we define S (S .... S,)T [ as the scaling vector. In scRNA-seq data, S; corresponds to the library size of the ith cell. The library size is related to the total sequencing reads and can be estimated by the sum of counts within each cell other methods 12 1 21. The (Xg(Xpp)
log(X, X,) = (oog(X.)....
follows multivariate normal random vector with mean 14"
and covariance Σr The precision matrix O = (E-))² represents the network. In summary, we have the following graphical model for count data, for 1 < i n.
Y; II Poisson (S S,),;,,
log X N (H', (0")-))
We denote this PLN distribution by Y. PLN(S; 11" E²). The above PLN model little different from the classical form 11. in which S; = 1 for 1 S n. We assume that the network O sparse and therefore we could maximize the penalized log- likelihoooo estimate O". However, the likelihood function in the PLN model involves a p-dimensional integration and is difficult to be computed, especially when p is large. Chiquet et al. developed a variational algorithm called VPLN to maximize the penalized log-likelihood [5] Although the variational method is computationally more feasible than directly maximizing the penalized likello,,,,,d, the estimator's theoretical properties difficult to be obtained.
We develop two-step estimator called PLNet that is computationally efficient and has good theoreti- cal properties. We first use the MMLE method to derive an estimator E for Σ* Then, we apply the D-trace method (30 with the covariance estimator E to estimate the sparse precision matrix O We show that the derived estimator θ is consistent estimator of O* even when the model dimension goes infinity exponentially as the sample size increases.
2.3 The maximum marginal likelihood estimator
We estimate μ = [uğ and E = [ojulesssss by maximizing the marginal log-likelinee log-likelihood to avoid high dimensional integrations. Let Y., be the jth column of Y. For S S P. we estimate nz and oj by maximizing the marginal log-likelihood of Y., and denote the estimator as ūy and وز and estimate ofx by maximizing the marginal log- likelihood of (Y.,, .y, Y.,) and
3

---

denote the estimator as 5k. The specific forms of these marginal log- li-liikiiioods and optimization problems are shown in the Supplemental Material (S2.1). Further, we define the MMLE for μ²" Σ as ji = [ Σ ( We apply the Newton- Paatonn-apphonnn algorithm to maximize the marginal log-likelihood with the initial values taken as the moment estimator fTT and nnTH in the Supplemental Material (S1.1). The explicit forms of the first and second-order partial derivatives of the marginal log- ligelike function are shown in the Supplemental Material (S2.1). The integrations involved in optimizations can be well approximated by the method described in [1]
The above MMLE Σ may be not positive semi-definite. However, the D-trace method requires the input covariance matrix estimator to be positive semi-definite to guarantee the convexity of the loss function. We project E to the space of positive semi-definite matrices and identify Σ that is closest to Σ in the space, i.e.,
= arg min A=
A
Σ
(2)
Then we add small positive definite matrix Σ and get positive-definite estimator E of Σ*
3 - ヒ || Σ
where I is the identity matrix. The optimization problem 2 for Σ can be solved by splitting conic solver [8]
With the covariance matrix estimator E. we apply the D-trace method to estimate the precision matrix,
1 (Io²) - tr (0) + X,, ||0ll,,or 2
arg min o-
(4)
The above optimization problem can be solved by an alternating direction method of multipliers
30 In this paper, we use a more efficient algorithm developed in [22] to calculate θ. We can also incorporate the prior information about the gene regulatory relationship to improve the performance of network inference. Let E be the set of gene pairs that cannot have direct interactions from the prior information, we can incorporate the prior information E by considering the following constrained optimization problem. o
{ tr (De?) - tr (9) + An ||0|1, 1,oft Oi = 0 for (i.j) E } 2
= arg min ez D
The tuning parameter An is selected by minimizing the following approximate Bayesian information criterion (BIC) (3I.
|655+ + 5 1 + (1111, log /n.
(5)
2
3 Theoretical properties
In this section, we establish the theoretical properties in the high dimensional setting. We first prove that E is a consistent estimator of Σ*. The convergence rate for E is similar to that of the sample covariance matrix for random variables with sub-Gaussian distribution. Then, under the same inrepresentability condition of the D-trace method in [30]. we derive the edge recovery property and consistency for the PLNet estimator θ.
3.1 Notation
Let G ((i,))) # 0} be positions of non-zero elements in O", G be the complement set of G, be the maximum node degree O*, 8 = ||0* and 0min = min |0;] be the minimal matimal absolute value of nonzero elements of O". Let Amax (A) and Amin (A) be the largest and smallest
eigenvalues of a symmetric matrix A. We define I r(E*) (E+ @ + @ E*) /2. )
where A 2 B is the Kronecker product. For a p X p matrix A, the row (i 1)p + j and column (k- - 1)p + 1 of (A) is I' (A) = (Ax1()=)) + Axx((i=))) /2. For two subsets Ti and T2

---

of {(i,)) $i.j S p), we define I (AA), Tu be the submatrix of I (A) whose rows and columns indexed by Ti and T2, respectively. Other notations are as follows,
= 1 -
max (KJ)EGr
kr (TT2,a) 1,00 ky Σ | 1,,00
3.2 The irrepresentability condition and rate of convergence
We first present the necessary irrepresentability condition for establishing the rate of convergence for the estimator in PLNet. This irrepresentabiiity condition is from the D-trace method in [30].
Condition 1 (Irrepresentability condition) There exists 0 < < 1, such that
11, - -a
max (KJ)EG
The irrepresentability condition I is equivalent to 2 a > 0. We also need a boundedness condition in the PLN model
Condition 2 (Boundedness condition) There exist positive constants M1, M2, My. M4. and Ms. such that the distribution of S; has bounded support [M1, M2]. maxigiAsp | G5 } M3
and M4 & Amin (E*) S Amar S Ms-
Based on the boundedness condition 2 we can establish the convergence rate for the covariance matrix estimator Σ.
Theorem (Rate of convergence for Σ) Under Condition 2 there exist positive constants A, B, such that for any
pr (15- E' 2() S p² Aexp(- BBne).
After plugging E. into the lasso penalized D-trace loss, get consistent estimator θ that converges to O* in several matrix norms.
Theorem (Rate of convergence for θ) Under the irrepresentability condition and the bounded- ness condition there exist constants A, B, such that for some 1) if
{12 (kg + kF)) + 5dkf) m B B (nlop log A max 12dkr, 12y (kgk? kr),
min (ss//,dd11 d (12 (ks:: + kG)) + 5dk?} > (0r) ( 1²
and
An 12 (KEk? kr) B
log n)
then with probability
22--
6 0r
5 (127
(kr + kf) + 5dk:3) B 1/2 logp + log A) 1/2 n 1/2
@ 0+
< 1/2 (127 (kisshi + 12)) + 5ck B 1/2 112 log log 1) 1/2
6 0*
min { 8 1/2 d + 1 (12 (Ke x 5dk B 1/2 (1 log log 1)2 1/2
Meanwhile, with a high probability, the sign of the sparse precision matrix O* can be recovered by 6. We have the following theorem about the sign consistency of é.
5

---

Theorem (Sign consistency for ①) Under the conditions in Theorem 2 for some 1j > 2, choosing the same and Xn in Theorem then with probability 1 pp²-", θ recovers all zeros and nonzeros in Or.
The rate-of-convergence and sign consistency results in Theorem 2 and 3 are closely related to the
tail probability of the MMLE estimator Σ in Theorem The boundedness condition is assumed to guarantee the convergence of the MMLE. Ignoring the complicated constants theorems, for any > 2, if we have n > Cnj log p. or in other words, if tends to infinity not faster than exponential of the sample size n. θ is consistent estimator of O Especially, the rate of convergence for θ is O (In(logp)/n) under Ise -nomm
4 Simulation studies
4.1 Simulation settings
We conduct simulations evaluate the performance of PLNet and compare it with the available network inference methods including VPLN [. glasso 77. and an estimator called PLNet-MOM by plugging the moment estimator 5nh into D-trace loss (See Supplemental Material (S1.1)). Both PLNet and VPLN are designed to estimate the precision matrix for count data in the PLN model. The glasso algorithm is classical approach for continuous data in GGM and apply glasso to the logarithmic transformation of the normalized data, which is defined Y,, (Y,, 1)/S,, 1 kis S n, where S; is the estimated library size of ith cell. all simulations, we estimate the library size by the total sum scaling (i.e. S = D j=1 Yij), which is a classical normalization method for scRNA-seq and is defined as the sum counts within each cell.
We consider 48 different simulation scenarios, which are 2 sample size setups (n = 500, 2000) x 3 dimension setups (p- = 1 00, 300, 500) x dropout levels (low and high) x 4 graph structures and simulate count data from the PLN model, where the dropout level represents the proportion of zeros the count matrix. Details of the simulation settings are shown the Supplementary Material (S3.1). In each scenario, independently repeat simulations 100 times. The four graph structures are as follows.
1. Banded Graph: Pairs (i.j) of nodes are connected if | jl S 2, i A All nonzero edges are set as 0.3.
2. Random Graph: Pairs of nodes are connected with probability 0.1. The nonzero edges set as 0.3 with probability 0.8 and as -0.3 with probability 0.2.
3. Scale-free Graph:
The Barabasi- Albert model 3 is used to generate a scale-free graph with
power 1. The nonzero edges are set as 0.3.
4. Blocked Graph:
The nodes are divided into blocks of equal sizes. Pairs of nodes in the same block are connected with probability 0.1 and the nonzero edges are set as 0.3. Different blocks have no edge connection.
The diagonal elements of precision matrices are set as 1 plus a positive number to guarantee positive definiteness.
4.2 Performance comparison
Table I and Table 2 show the area under the precision and recall curve (AUPR) of each estimator. AUPRS are calculated by varying the tuning parameters (i.e. the penalty parameters of the four estimators). As expected, the AUPR decreases as the number of genes increases. AUPRs in the high-dropout cases are generally smaller than that in the low-dropout cases. PLNet is the most robust estimation among these four estimators and outperforms PLNet-MOM. VPLN, and glasso in most simulation settings in AUPR, especially for the settings with high dropouts. For example, for 2000. p 100, Petet achieves an AUPR of 0.88 for the random graph under the scenario of the high dropout rate, while MOM, VPIN, and glasso only have AUPRs of 0.82, 0.69, and 0.18, respectively. Furthermore, PLNet-MOM outperforms glasso most settings and outperforms

---

Table I: Comparisons of PLNet, VPLN, glasso, and PLNet-MOM in terms of the area under precision and recall curve (AUPR) on simulation results for 500. The results are averages over 100 replicates with standard deviations in brackets.
Sample size
500
500
Ye 500
Dimension Dropout
p 100
p 300
p 500
Low
High
Low Banded graph
High
Low
High
PLNet
0.79 (0.03) 0,51 (0.04)
0.6 (0.06) 0.31 (0.04)
0,44 (0.08) 0.22 (0.03)
PLNet-MOM
0.62 (0.03) 0.43 (0.03)
0.33 (0.02) 0.21 (0.01)
0.21 (0.02) 0.13 (0.01)
VPIN
0.68 (0.06) 0.44 (0.04)
0.6 (0.02) 0.27 (0.02)
0.54 (0.01) 0.21 (0.02)
glasso
0.34 (0.04) 0.05 (0.01)
0.44 (0.03) 0.04 (0.01)
0.44 (0.02) 0.04 (0.01)
Random graph
PLNet
0.73 (0.07) 0.49 (0.07)
0.6 (0.07) 0.25 (0.05)
0.53 (0.09) 0.16 (0.04)
PLNet-MOM
0.6 (0.07) 0.44 (0.07)
0.44 (0.04) 0.22 (0.04)
0.36 (0.04) 0.14 (0.03)
VPIN
0.59 (0.07) 0.44 (0.08)
0.53 (0.09) 0.26 (0.05)
0.51 (0.06) 0.18 (0.05)
glasso
0.41 (0.05) 0.14 (0.02)
0.48 (0.05) 0.11 (0.02)
0.48 (0.05) 0.09 (0.02)
Scalc-free Graph
PLNet
0.73 (0.12) 0.55 (0.08)
0.67 (0.04) 0.44 (0.05)
0.59 (0.04) 0.35 (0.06)
PLNeI-MOM
0.63 (0.08) 0.51 (0.06)
0.48 (0.03) 0.34 (0.03)
0.39 (0.01) 0.27 (0.03)
VPLN
0.67 (0.05) 0.53 (0.09)
0.61 (0.03) 0.43 (0.04)
0.57 (0.07) 0.38 (0.05)
glasso
0.56 (0.06) 0.34 (0.05)
0.59 (0.02) 0.32 (0.03)
0,58 (0.02) 0.3 (0.04)
Blocked graph
PLNet
0.68 (0.05) 0.47 (0.08)
0.59 (0.08) 0.24 (0.07)
0.45 (0.1) 0.17 (0.04)
PLNet-MOM
0.58 (0.04) 0.43 (0.08)
0.43 (0.06) 0.2 (0.05)
0.32 (0.05) 0.15 (0.03)
VPIN
0.59 (0.04) 0.45 (0.07)
0.52 (0.07) 0.24 (0.06)
0.44 (0.08) 0.17 (0.04)
glasso
0.37 (0.03) 0.16 (0.03)
0.44 (0.06) 0.11 (0.03)
0.4 (0.06) 0.09 (0.02)
Table 2: Comparisons of PLNet, VPLN, glasso, and PLNet-MOM in terms of the area under precision and recall curve (AUPR) on simulation results for n 2000. The results are averages over 100 replicates with standard deviations in brackets.
Sample size
n = 2000
2000
2000
Dimension
p 100
p 300
p 500
Dropout
Low
High
Low High
Low
High
Banded graph
PLNct
0.99 (0.01) 0.96 (0.01)
0.99 (0.01) 0.94 (0.02)
0.98 (0.01) 0.89 (0.08)
PLNet-MOM
0.97 (0.01) 0.92 (0.01)
0.91 (0.01) 0.83 (0.02)
0.83 (0.02) 0.75 (0.01)
VPLN
0.95 (0.01) 0.89 (0.03)
0.94 (0.01) 0.79 (0.15)
0.94 (0.01) 0.81 (0.01)
glasso
0.62 (0.03) 0.04 (0.01)
0.82 (0.01) 0.07 (0.01)
0.85 (0.01) 0.15 (0.02)
Random graph
PLNet
0.98 (0.01) 0.88 (0.04)
0.98 (0.03) 0.85 (0.05)
0.99 (0.90) 0.83 (0.04)
PLNet-MOM
0.94 (0.02) 0.82 (0.06)
0.94 (0.01) 0.77 (0.05)
0.93 (0.01) 0.74 (0.05)
VPLN
0.78 (0.08) 0.69 (0.07)
0.88 (0.03) 0.67 (0.1)
0.86 (0.11) 0.67 (0.11)
glasso
0.55 (0.06) 0.18 (0.03)
0.8 (0.03) 0.24 (0.04)
0.84 (0.02) 0.26 (0.04)
Scale-free Graph
PLNet
0.89 (0.17) 0.85 (0.11)
0.97 (0.02) 0.85 (0.03)
0.96 (0.03) 0.83 (0.02)
PLNet-MOM
0.85 (0.11) 0.81 (0.08)
0.86 (0.01) 0.75 (0.02)
0.83 (0.01) 0.71 (0.01)
VPLN
0.74 (0.16) 0.67 (0.15)
0.79 (0.04) 0.68 (0.11)
0.8 (0.05) 0.66 (0.13)
glasso
0.59 (0.14) 0.45 (0.06)
0.78 (0.02) 0.5 (0.03) Blocked graph
0.81 (0.02) 0.53 (0.02)
PLNet
0.94 (0.02) 0.83 (0.07)
0.97 (0.01) 0.81 (0.08)
0.97 (0.01) 0.77 (0.05)
PLNet-MOM
0.88 (0.04) 0.75 (0.08)
0.91 (0.02) 0.72 (0.08)
0.89 (0.02) 0.68 (0.05)
VPLN
0.73 (0.03) 0.66 (0.07)
0.78 (0.04) 0.62 (0.07)
0.8 (0.06) 0.59 (0.11)
glasso
0.47 (0.05) 0.2 (0.03)
0.7 (0.04) 0.21 (0.04)
0.75 (0.03) 0.21 (0.04)
7

---

VPLN under the scenario of large sample sizes. The results of the area under the receiver operating characteristic curve (AUC) are similar and shown in Supplementary Table S1-S2.
To further demonstrate the performance of PLNet, we visualize the mean networks predicted by the four methods tuned by BIC criterion for the banded graph with n 2000, p 100 over the 100 simulations (Figure []. More specifically,, we calculate the relative frequency Fij that an algorithm reports edges over the 100 simulations. For positions (i,j) with O. * Púj is the proportion that an algorithm correctly recovers the edge in the 100 simulations, while for positions (i. j) with O2 F.y is the proportion that an algorithm falsely predicts edges between nodes i and j in the 100 simulations. We plot the relative frequency matrices of the four methods in Figure [ We clearly see that PLNet is able to detect more true positives while having fewer false positives than other methods, especially for the high dropout case. The results of mean predicted networks for the other three graphs with n 2000, p 100 are similar and shown Supplementary Figure $1-S3.
PLNct
PLNct MOM
VPLN
glasso
Inodolp MD1
Truth
Frequancy
-0.5
poodog 4SH
-1.0
Figure 1:
The mean networks predicted by PLNet, VPLN, glasso, and PLNet-MOM for the banded graph with 100 nodes and n = 2000. False edges are colored in red and true edges are in blue. The left panel is the network matrix for reference.
Supplementary Table S3-S4 show the computational time of the four methods. PLNet is computa- tionally more efficient than VPLN in most simulation settings. Also, the computational efficiency of PLNet is roughly the same for different sample sizes, while the computational complexity of VPLN maintains a linear relationship with sample size. The detailed comparisons between the computational efficiency of the four methods are shown in the Supplementary Material (S3.2.2).
5 Application to scRNA-seq dataset
We apply PLNet and VPLN to infer the gene regulatory network of CD14+ Monocytes profiled in Kang al. [TO] The single cells are profiled in two different conditions, IFN-3-treetedd and controlled. IFN-3 is a cytokine in the interferon family that influences the transcriptional profiles for many genes, especially those in the JAK/STAT pathway [16]. We focus on the IFN-fttreated cells (2147 cells) and use the top 200 highly variable genes that are used in (20 for network analysis.
We compare networks of PLNet and VPLN with parameters tuned such that the network densities are around 5%. Gene Ontology (GO) analysis [1] shows that the 200 genes mainly involve in 4 biological processes, including immunity" (Module M2). "Cellular protein metabolic process" (Module My), and "Proteolysis" (Module M4). FigureZ the predicted networks of the genes in the 4 modules by Pl.Net and VPLN, where the colors represent the partial correlations between genes. The partial correlation
mediated signaling pathway" (Module MD). "Neutrophil-mmdiitte
given by PLNet between genes i and j is defined as - given by VPLN is defined similarly. We clearly see that the network given by PLNet tends to have more connections within the modules than VPLN. To see this more clearly, for each module Mr. we calculate the ratio between within-module and between-module connections R(Mss) = Ei Wiz/ W ExEM,, jeMe Wij. where the weights Wiy are set as the absolute partial correlation between genes i and j (weighted within-between connection ratio) or are set as and 0 depending on whether genes and j are connected (unweighted within-between connection ratio). The within-
The partial correlation
8

---

(a) PLNet
(b) VPLN
Figure 2: Heat maps of partial correlations between genes in the 4 GO modules given by PLNet (a) and VPLN (b). Red: Cytokine-mediated signaling pathway (Module M,); Orange: Neutrophil- mediated immunity (Module M2); Green: Cellular protein metabolic process (Module M3); Blue: Proteolysis (Module M4).
between connection ratios PLNet are much larger than that of VPLN most cases (Table 3) Similar results also hold for networks of other different densities that are shown in the Supplementary Material (S4.1).
Table The within- betteen connection ratios of the modules in the networks estimated by PLNet and VPLN tuned such that the network densities are around 5%
Type Weighted
Method
Module M1 Module M2 Module M3 Module M4
PLNet
0.681
0.231
0.582
0.387
VPLN 0.621
0.381
0.459
0.259
Unweighted
PLNe 0.533
0.129
0.419
0,383
VPLN 0,515
0.200
0.170
0.213
To further compare the performance of PLNet and VPLN the real scRNA-seq data analysis, we construct silver standard based on available regulatory network database obtained from ChIP-seq experiments (the hTFtarget database (29 The silver standard consists of edges between the 200 highly variable genes and 36 transcription factors (TF) from the top 500 highly variable genes. We consider additional TFs from top 500 highly variable genes because the top 200 highly variable genes only contain 10 TFs and the silver standard constructed only using these genes is too sparse to a reliable silver standard. We compare networks estimated by PLNet and VPLN with their parameters tuned such that the network densities at the same levels (1% to 10%). Table 4 shows the number true edges in the silver standard detected by PLNet and VPLN at different network densities. Clearly, see that PLNet detects more true edges at the same density levels than VPLN, suggesting that PLNet has a higher true discovery rate than VPLN.
Table The number true edges estimated by two methods with different density levels.
Density 0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
PLNet
16
23
35
41
44
62
73
81
92
VPLN
5
20
48
Discussion
In this paper, we consider the PI.N graphical model for count data. This model has an intuitive explanation for single-cell gene regulatory network analysis. The Poisson layer is capture the

---

technical noises and the log- log-normal layer is to model the biological fluctuations in single cells. Gene regulatory network is represented by the precision matrix of the latent log-normal model. To estimate the precision matrix, we propose two-step estimator PLNet, using the MMLE to estimate the covariance matrix and then minimizing the penalized D-trace loss to estimate the precision matrix. The simplicity of this estimation procedure allows us to establish the consistency theory for the proposed PLNet estimator for the high dimensional setting. The numerical analysis also shows that the PLNet method outperforms available methods.
The proposed method can be generalized in several ways. First, we may consider the zero-inflated PLN (ZIPLN) model to account for the excessed number of zeros that cannot be modeled by PLN models. We could also use the MMLE method to estimate the covariance matrix of the ZIPLN model, but the convergence theory is more difficult to develop. Secondly, straightforward generalization the differential network analysis. We could use the same method to estimate covariance matrices in two different conditions and use the D-trace loss developed in earlier work [27] for differential network analysis in single cells. Thirdly, the mean parameters of the latent log-normal random variable can depend on covariates. In such model, we could first use regression to estimate the mean parameters and derive the corresponding moment estimator of the covariance matrix. Another generalization is gene regulatory network analysis of mixtures of cell populations. Different cell populations may have different gene regulatory networks and we could jointly model the mixture and infer the gene regulatory networks for all cell populations. Finally, the current method is based on parametric inference. If the model is mis-specified network inference could be misleading. To mitigate this problem, normal model with a Gaussian copula graphical model, but how to do network inference with this semi-parametric model is still challenge that needs to be overcome in the future.
consider a semi parametrri model by replacing the latent log-
Acknowledgments and Disclosure of Funding
We thank Changhu Wang and Siyuan Huang for helpful discussions. This work was supported by the National Natural Science Foundation of China [No. 11971039 to R.X.J. the National Key Basic Research Program of China [2020YFE0204000 to R.X.J, and Sino-Russian Mathematics Center.
References
[1] Aitchison and C. H. Ho. The multivariate Poisson-log normal distribution. Biometrika, 76(4): 643-653, 1989.
[2] Genevera I. Allen and Zhandong Liu. A local Poisson graphical model for inferring networks from sequencing data. IEEE Transactions on NanoBioscience, 12(3):189-198, 2013.
[3] Albert-Lász1ó Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286(5439):509-555 1999.
[4] Martin Bengtsson, Anders Ståhlberg, Patrik Rorsman, and Mikael Kubista. Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels. Genome Research, 15(10):1388-133 2005.
[5] Julien Chiquet, Stephane Robin, and Mahendra Mariadassou. Variational inference for sparse network reconstruction from count data. In International Conference on Machine Learning. pages 1162-1171. PMLR, 2019.
[6] Jeffrey A. Farrell, Yiqun Wang, Samantha J. Riesenfeld, Karthik Shekhar, Aviv Regev, and Alexander Schier. Single-cell reconstruction of developmental trajectories during zebrafish embryogenesis. Science,
2018.
[7] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432-441, 2008
[8] Anqi Fu, Balasubramanian Narasimhan, and Stephen Boyd. CVXR: An R package for dis- ciplined convex optimization. Journal of Statistical Software, 94(14):1-34, 2020. ISSN 1548-7660.
10

---

[9] David I. Inouye, Eunho Yang. Genevera I. Allen, and Pradeep Ravikumar, A review of multi- variate distributions for count data derived from the Poisson distribution. Wiley Interdisciplinary Reviews:
Computational Statistics, ((3)=1138, 2017.
[10] Hyun Min Kang, Meena Subramaniam, Sasha Targ. Michelle Nguyen, Lenka Maliskova, Elizabeth McCarthy, Eunice Wan, Simon Wong, Lauren Byrnes, Cristina M. Lanata, Rachel E. Gate, Sara Mostafavi, Alexander Marson, Noah Zaitlen, Lindsey A. Criswell, and Chun Jimmie Ye. Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology, 36(1):89-94, 2018.
[11] Maxim V Kuleshov, Matthew Jones, Andrew D Rouillard, Nicolas F Fernandez, Qiaonan Duan, Zichen Wang, Simon Koplev, Sherry I. Jenkins, Kathleen M Jagodnik, Alexander Lachmann, et al. Enrichr: comprchensive gene set enrichment analysis web server 2016
update. Nucleic Acids Research, 44(W1): 2016.
[12] Michael I Love, Wolfgang Huber, and Simon Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12):1-21, 2014.
[13] Aaron TL Lun, Karsten Bach, and John C Marioni. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biology, 17(1):1-14, 2016.
[14] Shisong Ma. Qingqiu Gong, and Hans Bohnert. An Arabidopsis gene network based on the graphical Gaussian model. Genome Research, 17(11):1614-1625, 2007.
[15] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable selection with the Lasso. The Annals of Statistics, 34(3):1436-146) 2006.
[16] Sara Mostafavi, Hideyuki Yoshida, Devapregasan Moodley, Hugo LeBoité, Katherine Rothamel, Towfique Raj. Chun Jimmie Ye, Nicolas Chevrier, Shen-Ying Zhang, Ting Feng, Mark Lee, Jean-Laurent Casanova, James D. Clark, Martin Hegen, Jean-Baptiste Telliez, Nir Hacohen, Philip L. De Jager, Aviv Regev, Diane Mathis, and Christophe Benoist. Parsing the interferon transcriptional network and its disease associations. Cell, 164(3):564 578, 2016.
[17] C. Naylor and A. M. Smith. Applications of method for the efficient computation of posterior distributions. Journal of the Royal Statistical Society. Series C (Applied Statistics), 31 (3):214-22) 225, 1982.
[18] Jie Peng, Pei Wang, Nengfeng Zhou, and Ji Zhu. Partial correlation estimation by joint sparse regression models. Journal of the American Statistical Association, 104(486):735-746, 2009.
[19] Mark D. Robinson, Davis McCarthy, and Gordon K. Smyth. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1): 139-140, 2009.
[20] Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister, Efthymia Papalexi, William M. Mauck, Yuhan Hao, Marlon Stoeckius, Peter Smibert, and Rahul Satija. Compre- hensive integration of single-cell data. Cell,
2019.
[21] Catalina A Vallejos, John C Marioni, and Sylvia Richardson. BASiCS: Bayesian analysis of single-cell sequencing data. PLOS Computational Biology, 11( 1 004333, 2015.
[22] Cheng Wang and Binyan Jiang. An efficient ADMM algorithm for high dimensional precision matrix estimation via penalized quadratic loss. Computational Statistics Data Analysis, 142: 106812, 2020.
[23] Anja Wille, Philip Zimmermann, Eva Vranová, Andreas Fürholz, Oliver Laule, Stefan Bleuler, Lars Hennig, Amela Prelić, Peter von Rohr, Lothar Thiele, et al. Sparse graphical Gaussian modeling of the isoprenoid gene network in Arabidopsis thaliana. Genome Biology, 5(11):1-13, 2004.
[24] Hao Wu, Xinwei Deng, and Naren Ramakrishnan. Sparse estimation of multivariate Poisson log-normal models from count data. Statistical Analysis and Data Mining: The ASA Data Science Journal, 11(2):66-77, 2018.

---

[25] Eunho Yang, Pradeep K Ravikumar, Genevera I Allen, and Zhandong Liu. On poisson graphical models. In Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
[26] Jianxin Yin and Hongzhe Li. A sparse conditional Gaussian graphical model for analysis of genetical genomics data. The Annals of Applied Statistics, 5(4):2630, 2011.
[27] Huili Yuan, Ruibin Xi, Chong Chen, and Minghua Deng. Differential network analysis via lasso penalized D-trace loss. Biometrika, 104(4):755-770, 2017.
[28] Ming Yuan and Yi Lin. Model selection and estimation in the Gaussian graphical model. Biometrika, 94(1):19-35, 2007.
[29] Qiong Zhang, Wei Liu, Hong-Mei Zhang, Gui-Yan Xie, Ya-Ru Miao, Mengxuan Xia, and An- Yuan Guo. htftarget: comprehensive database for regulations of human transcription factors and their targets. Genomics, Proteomics & Bioinformatics, 18(2):120-128, 2020.
[30] Teng Zhang and Hui Zou. Sparse precision matrix estimation via lasso penalized D-trace loss. Biometrika, 101(1):103-125) 2014.
[31] Sihai Dave Zhao, T. Tony Cai, and Hongzhe Li. Direct estimation of differential networks. Biometrika, 101(2):253-265 2014.
[32] Grace XY Zheng, Jessica M Terry, Phillip Belgrader, Paul Ryvkin, Zachary W Bent, Ryan Wilson, Solongo B Ziraldo, Tobias D Wheeler, Geoff P McDermott, Junjie Zhu, et al. Massively parallel digital transcriptional profiling of single cells. Nature Communications, 8(1):1-12. 2017.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Section 6
(c) Did you discuss any potential negative societal impacts of your work? [No] Our work do not have any negative societal impact.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Conditions aadd Theorems.
(b) Did you include complete proofs of all theoretical results? [Yes] See Supplemental Material for proofs.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main exper- imental results (either in the supplemental material or as URL)? [Yes] See Supple- mental Material for codes.
(b) Did you specify all the training details (e.g.. data splits, hyperparameters, how they were chosen)? [Yes] See Section 4.1 for Simulation settings and equation for hyperparameter chosen.
(c) Did you report error bars (e.g., with respect to the random seed after running exper- iments multiple times)? [Yes] See the standard deviations Table and for error bars.
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Supplemental Material (S3.2.2).
4. If you using existing assets (c.g., code, data, models) or curating/releasnng assets...
12

---

(a) If your work uses existing assets, did you cite the creators? [Yes] See lines 105, 1 114,
128, 135, 138, 177, 195, 230 for citations.
(b) Did you mention the license of the assets? [Yes] All the assets are open available.
(c) Did you include any new assets either in the supplemental material or as URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [No]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [No]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals. if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "xwPWfbye_Pg", "reviewer": "Reviewer_gBop", "summary": "The authors propose a new algorithm (PLNet) for estimating the precision matrix of a Poisson log-normal (PLN) graphical model for count data. The authors evaluate their method on both simulated and real-world data and find improvements over baseline methods.", "questions": "1. Before reviewing this manuscript, I was not familiar with the VPLN model. As such, it wasn't clear to me why PLNet is computationally more efficient than VPLN. It would be great if the authors could provide more details on this.", "limitations": "The authors addressed the limitations of their work in a satisfactory manner. They did not discuss potential negative societal impact, which is fine for this work (I do not see any potential negative impacts here).", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall I was quite happy with the paper and recommend acceptance.\n\nStrengths:\n\n1. The authors tackle a high-impact problem; rich scRNA-seq datasets are becoming more and more common, and methods for analyzing them are in high demand. Moreover, many works have found that analyzing scRNA-seq datasets after first normalizing and log-transforming the data leads to significant distortions. As such, methods that are able to operate on raw counts are crucial.\n2. To my knowledge, the authors' proposed algorithm (as well as the theoretical analysis of it) is novel.\n3. The authors demonstrate that their algorithm outperforms naive GGM approaches as well as a PLN approach that optimizes a lower bound to the log-likelihood (as opposed to directly optimizing the likelihood).\n4. The writing in the manuscript was clear and I generally found the paper to be well-structured.\n\nWeaknesses:\n\n1. Ideally I would've liked to see another application on real-world data, though given space constraints I understand that may not have been possible to fit in.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Yg9HRXQ4idG", "reviewer": "Reviewer_zkNW", "summary": "This paper proposes PLNet, a two-step approach to estimate the precision matrix of the latent normal distribution which represents the graph of count data. PLNet first considers the Poisson log-normal (PLN) graphical model and utilizes the maximum marginal likelihood estimator (MMLE) to estimate the latent covariance matrix. In the second step, PLNet estimates the sparse precision matrix by minimizing the lasso penalized D-trace loss. The authors also prove that the resulting estimator is a consistent estimator of the precision matrix in high dimensional settings under mild conditions. Extensive experiments, especially simulation experiments, demonstrate that PLNet outperforms baseline methods in both accuracy and efficiency.", "questions": "1. The reviewer wonders if it is possible to incorporate the gene regulatory network information in the modeling process to improve the estimation performance?\n\n2. In scRNA-seq data, sometimes gene expression of some genes is not detected. The reviewer wonders if the authors consider using a probability describing whether a gene's expression is detected to better approximate the real data.", "limitations": "Authors have discussed the Limitations and societal impacts in the Conclusion Section.\n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. Instead of using Gaussian distribution, the authors leverage multivariate Poisson log-normal (PLN) distribution to estimate the graphical models of count data, which handles the over-dispersion problem and is more likely to approximate the distribution of count data like scRNA-seq data. \n\n2. PLNet is computationally efficient as it avoids computing high-dimensional numeric integration when computing the log-likelihood of the multivariate PLN model, and applying the D-trace method to estimate the precision matrix is more efficient.\n\n3. The authors provide the theoretical analysis of the estimator in the high-dimensional setting and prove that the estimator is a consistent estimator of the precision matrix under mild conditions.\n\nWeaknesses:\n\n1. It seems that the authors did not mention how to obtain the mean vector of the latent normal distribution in the PLNet. It would be great if the authors could discuss it.\n\n2. Sometimes in scRNA-seq data, the number of genes (i.e., $p$) may be large, and thus the reviewer wonders if the estimator is still a consistent estimator of the precision matrix in this case?\n\n3. Although extensive simulation results have been presented in this paper, it seems that the real-world experimental results are not very clear. The reviewer thinks that the authors may use some other criteria to quantify the results rather than just using the heatmap and within-between connection ratios. \n\n4. The presentation of experiment results should be improved, and it might be better to use tables to present the results.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "IU2MvaT4WI", "reviewer": "Reviewer_CQBu", "summary": "This paper introduces a consistent estimator for the sparse precision matrix of the Poisson log-normal (PLN) graphical model. The paper provides theoretical guarantees of the estimator's consistency and convergence rate. Using simulated data, the paper demonstrates that the proposed estimator is more accurate than competitive baselines across a variety of conditions. The paper also compares the proposed estimator to a recent baseline on inferring a gene regulatory network from RNA sequence data.", "questions": "* In the experiments section there are 3 instances where PLNet and VPLN seem to achieve the same AUPR, but only PLNet is bolded. Is this just due to rounding?\n\n* What is the relationship of the proposed method to the method given in section 2.4.2 of Fan, Liao, and Liu (2015) [1]?\n\n[1] An Overview on the Estimation of Large Covariance and Precision Matrices. Jianqing Fan, Yuan Liao and Han Liu. 2015.\nhttps://arxiv.org/pdf/1504.02995.pdf", "limitations": "Yes.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Clarity\nThis paper is very well-written and clear. It provides excellent coverage of prior work and situates the proposed approach well within the relevant literature. The technical exposition is cleanly structured and easy to follow. The experiments are reported with sufficient detail. I also appreciate that the proposed method is grounded in an applied problem, which provides strong motivation for the paper, and helps in understanding the method.\n\n## Quality \nI found no flaws in the paper's logic or in the proofs for theorems 1-3 (which I only skimmed). The proposed method is simple and practical, while the theory guarantees root-n convergence of the estimator under mild conditions.\n\nThe experiments are thorough. The paper does a good job of simulating many synthetic data sets under varying conditions and performance to reasonable baselines. One drawback is that only the simple case of a banded covariance matrix is given as an example to compare visually the results of the proposed method to the baselines. Perhaps the authors can add more visual inspections to the appendix.\n\nMinor comments/edits:\n- \"semi-positive definite\" -> \"positive semi-definite\"\n- line 69: \"computational\" -> \"compuationally\"\n- line 104: \"ditribution\" -> \"distribution\"\n- line 170 \"convergency\" -> \"convergence\"\n- line 215 \"tunned\" -> \"tuned\"\n\n## Originality\nI've not worked in the area of Gaussian graphical models, nor in sparse estimation of precision/covariance matrices, so it's difficult for me to judge how original or novel the paper is. My sense is that the method itself is  built from standard pieces and is not very original, but the theory that establishes the method's consistency and convergence rates is non-trivial and new. \n\n## Significance\nWhile this paper focuses on the particular applied setting of inferring gene networks from RNA sequence data, the method is fairly general to any setting in which multivariate count data arises. The method is simple, easy to understand, and seems easy to implement, so I expect it will be popular among practitioners. I would not be surprised if similar methods are already in use, and thus the new theoretic guarantees that this paper contributes will likely make the proposed method popular in practice.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept