 Methodology article
 Open Access
 Published:
Network analysis for count data with excess zeros
BMC Genetics volume 18, Article number: 93 (2017)
Abstract
Background
Undirected graphical models or Markov random fields have been a popular class of models for representing conditional dependence relationships between nodes. In particular, Markov networks help us to understand complex interactions between genes in biological processes of a cell. Local Poisson models seem to be promising in modeling positive as well as negative dependencies for count data. Furthermore, when zero counts are more frequent than are expected, excess zeros should be considered in the model.
Methods
We present a penalized Poisson graphical model for zero inflated count data and derive an expectationmaximization (EM) algorithm built on coordinate descent. Our method is shown to be effective through simulated and real data analysis.
Results
Results from the simulated data indicate that our method outperforms the local Poisson graphical model in the presence of excess zeros. In an application to a RNA sequencing data, we also investigate the gender effect by comparing the estimated networks according to different genders. Our method may help us in identifying biological pathways linked to sex hormone regulation and thus understanding underlying mechanisms of the gender differences.
Conclusions
We have presented a penalized version of zero inflated spatial Poisson regression and derive an efficient EM algorithm built on coordinate descent. We discuss possible improvements of our method as well as potential research directions associated with our findings from the RNA sequencing data.
Background
Graphical models help us to explore relationships between nodes in graphs. Undirected graphical models or Markov random fields have been a popular class of models for representing conditional dependence relationships between nodes. Examples include Gaussian graphical models for continuous data, Ising model for binary data, and multinomial graphical models. These Markov networks help us to understand complex interactions between genes in biological processes of a cell and have been well studied in bioinformatics. Examples of Markov networks in learning the network structure from microarray and next generation sequencing data include [1–4]. For more details on Markov network inference, see those and the references therein.
The main focus of this study is to infer the network structure for a count data. The autoPoisson model in [5] is a natural extension of univariate Poisson distribution. However it can model only negative dependencies, so that the conditional distributions define a unique joint distribution consistently. Yang et al. [6] propose variants of the autoPoisson model such as truncated, quadratic, and sublinear Poisson graphical models(PGM). However none of them provide a satisfactory answer to the question of how to specify a consistent joint graphical model for count data capturing both positive and negative dependencies. Allen and Liu [4] consider a local PGM (LPGM). The LPGM does not have a consistent joint graphical model, but it has the local Markov property and thus the zero coefficient of an edge weight between two nodes implies the conditional independence of the two nodes given the others. žitnik and Zupan [7] consider a latent factor Poisson model and [8] propose to learn conditional dependence structures for binary and Poisson data via marginal loss functions. Also a semiparametric Guassian copula, called the nonparanormal graphical model (NPGM), has been proposed [9].
In practice, zero counts are sometimes more frequent than are expected under a univariate Poisson distribution. In such cases, a zeroinflated Poisson (ZIP) distribution is often adopted. Applications of ZIP models include modeling of defects in quality control [10] and alcoholism and substance abuse in medicine [11]. Extensions of a ZIP model in different frameworks are wellstudied in the literature. Dobbie and Welsh [12] extend the two component approach in [13] for serially correlated count data exhibiting extra zeros. Monod [14] develops a zeroinflated spatial Poisson (ZISP) model. Buu et al. [11] study variable selection methods such as LASSO and onestep SCAD for ZIP regression models. For computation, a local linear approximation (LLA) is adopted. The LLA algorithm fails to converge particularly with small sample sizes because it requires fitting unpenalized ZIP regression models. Wang et al. [15] propose an expectation maximization (EM) algorithm [16] for a penalized ZIP regression model built on coordinate descent algorithms. The EM algorithm seems to have some advantages over the LLA algorithm in numerical convergence and tuning.
In this paper, we are interested in the construction of graphical models for count data, particularly, with excessive zeros. To this end, we propose a penalized version the ZISP model in [14] called zero inflated local Poisson graphical model (ZILPGM) and derive an EM algorithm built on coordinate descent as in [15]. We show the effectiveness of our method on simulated and real data. In an application to a RNA sequencing data, we investigate the gender effect by comparing the estimated networks according to different genders. It has been well noted that gender is one of the major contributors in the differentiation of gene expression profiles [17, 18] and various sexually dimorphic phenotypes, most of which result from hormonal differences [19]. It was reported that transcriptome study could be predicted to represent a different promising approach for the identification of biological pathways linked to sex hormone regulation and the analysis of associated gene regulatory networks [20]. However, the elucidation of underlying mechanisms of the gender differences is still an area of interest and intense investigation.
The paper is organized as follows. In “Methods” section, we propose a new graph learning method based on ZISP and provide an efficient EM type numerical algorithm. In “Results” section, we compare performances of our method with LPGM on simulated and real data sets. Some discussions and concluding remarks are given in “Conclusions” section.
Methods
In this section, we present our graph learning method based on a penalization of the ZISP in [14] and derive an efficient EM algorithm for its computation.
Zero inflated local Poisson graphical model
Let N denote the number of observations and p denote the number of variables or nodes. Denote \(\mathcal {G}=(V, E)\), where V={1,…,p} is the set of vertices or nodes and E is the set of edges. We use uppercase letters such as X and Z when we refer to random variables. Observations are written in lowercase. For example, x _{ i } denote ith observation of X. Vectors and matrices are represented by boldface and blackboard boldface letters, respectively. Define \(\mathbb {X}=(x_{ij})_{N \times p}\), where x _{ ij } is generated from two latent components with zero and Poisson states. Let z _{ ij } be a latent variable such that z _{ ij }=1 if x _{ ij } is from zero state and z _{ ij }=0 if x _{ ij } is from Poisson state. z _{ ij } follows a Bernoulli distribution with π _{ j }. Let I(·) denotes an indicator function. Then the ZISP model in [14] is defined by
where \(\mu _{j} = \exp \left (\beta _{j} + \sum _{k\neq j} \beta _{jk}x_{k} \right)\), β _{ j } is an intercept adjusting for X _{ j }, and β _{ jk } is the parameter accounting for the conditional relation between X _{ j } and X _{ k }.
Due to the zero inflation term in the conditional probability, the situation becomes more complicated in our case than in LPGM. Because the important part is the pairwise interaction term in the pairwiseonly dependency models, the situation is basically similar. In order to have a valid joint distribution, the coefficient for the interaction term β _{ jk } should be nonpositive. As in the LPGM, we do not solve the issue of negative parameters in the Poisson graphical model. Note that any existing approaches (e.g. in [6]) do not succeed in giving a satisfactory answer to the consistency issue. Rather, we focus not on the consistency issue but on the practical issue of estimating positive as well as negative dependencies as in LPGM.
In order to learn graph structures, we consider the minimization of the penalized pseudo loglikelihood of (1) in the general weighted LASSO form:
where \(\mu _{ij} = \exp \left (\beta _{j} + \sum _{k\neq j} \beta _{jk}x_{ik} \right)\), λ≥0 is the penalty parameter, and w _{ jk }≥0 is an appropriate weight. As in [4], we can select the tuning parameter using the stability selection criterion in [21]. More specifically, we select the optimal λ is selected from 30 equalspaced grid points in log scale on [ λ ^{max},λ ^{min}], where \(\lambda ^{\text {max}}=\max _{j \in \{1,\cdots,p\}}\max _{k\neq j} \frac {1}{N} \sum _{i=1}^{N} x_{ik}x_{ij}\) and λ ^{min}=λ ^{max}×10^{−4}. For each j, we fit poisson regression using glmnet. Then w _{ jk }=1 for covariates with nonzero coefficients. Otherwise, w _{ jk } is set to be sufficiently large value, e.g., 10^{5}. Note that the purpose of the penalization is to select spatial neighbors. If β _{ jk }=0, then X _{ j } and X _{ k } is declared to be conditionally independent of the other variables.
The penalized pseudo loglikelihood in (2) is separable with respect to the coordinate index. Hence minimizing (2) is equivalent to separately minimizing the p coordinate functions:
Details on the algorithm is discussed later in this section. Once we solve (3), we can estimate the graph structure from the estimated set of edges: \(\hat {E} = \{ (j,k) : \hat {\beta }_{jk} \neq 0~\text {or}~\hat {\beta }_{kj} \neq 0, j \neq k\}\). We devise an EM algorithm as in [15] to minimize (3).
Computational algorithm
Let \(\mathcal {O}_{j}=\{i:x_{ij} = 0\}\) and \(\mathcal {P}_{j}=\{i:x_{ij} \neq 0\}.\) The negative loglikelihood function in (2) is the sum of
for j=1,…,p. However, it is difficult to maximize this likelihood directly because the score function of \(\sum _{i \in \mathcal {O}_{j}} \log \Big ({\pi _{j}}/\left (1\pi _{j}\right) + {e^{\mu _{ij}}} \Big)\) cannot be simplified [14, 22].
Instead of a direct optimization of the likelihood function, we express the likelihood function as a mixture distribution by introducing a latent variable and derive an EM algorithm.
Define β _{−j }=(β _{0},(β _{ k })_{ k≠j }))^{T} and x _{ i,−j }=(1,(x _{ ik })_{ k≠j }))^{T}. The loglikelihood function with respect to complete data can be written as
The decomposed likelihood function in the above can be easily maximized via an EM algorithm alternating between the expectation of the complete data likelihood over the latent variable z _{ ij } and the maximization of the likelihood given z _{ ij }’s.
Define the responsibility of zero state for jth variable on ith observation at mth step as \(z_{ij}^{(m)}=\mathbb {E} \left (z_{ij}x_{ij},\boldsymbol {\beta }_{j}^{(m)}\right)\) and the probability of zero state at mth step as
Our EM algorithm alternates the following steps until convergence.

Estep: Estimate z _{ ij } by its conditional mean \(z_{ij}^{(m)}\) given data and parameters from the previous step.
$$z_{ij}^{(m)}= \left\{ \begin{array}{ll} \frac{\pi_{j}^{(m)}}{\pi_{j}^{(m)}+\left(1\pi_{j}^{(m)}\right)\exp\left(\mu_{ij}^{(m)}\right)} & \text{if } x_{ij}=0,\\ 0, & \text{if } x_{ij}=1,2,\ldots \end{array}\right. $$ 
Mstep : Estimate \(\boldsymbol {\beta }_{j}^{(m)}\).
Here we set the initial values for our EM iteration as \(\pi _{j}^{(0)}=\) the number of zeros of jth variable /n for j=1,⋯,p and \(\boldsymbol {\beta }_{j}^{(0)}=\mathbf {0}\).
Now let us discuss the estimation of \(\boldsymbol {\beta }_{j}^{(m)}\) in detail. For each variable, we use the MajorizeMinimization (MM) algorithm in [23], which extends the central idea of EM algorithms to situations not necessarily involving missing data nor even maximum likelihood estimation. A function g(θθ _{ m }) is said to majorize a function f(θ) at θ _{ m } provided that f(θ _{ m })=g(θ _{ m }θ _{ m }) and f(θ)≤g(θθ _{ m }) for θ≠θ _{ m }. The key idea is that the surrogate majorizing function g(θθ _{ m }) is minimized iteratively, instead of the original objective function f(θ) with the nonquadratic log likelihood and the nondifferentiable sparsity inducing penalty [23]. The MM algorithm starts from an initial guess, θ _{0}. Let θ _{ m+1} denote the minimizer of the surrogate g(θθ _{ m }). Then the following inequalities hold:
The above inequality can easily be shown by definition of θ _{ m+1} and the majorization conditions. The descent property makes the MM algorithm numerically stable [24].
The objective to maximize is \(l_{j}^{c2} = \sum _{i=1}^{N}\left (1z_{ij}\right) \left (x_{ij} \boldsymbol {x}_{i,j}^{T}\boldsymbol {\beta }_{j}\exp \left (\boldsymbol {x}_{i,j}^{T}\boldsymbol {\beta }_{j}\right)\log x_{ij}!\right)\) whose first and second derivatives with respect to β _{−j } are
Let X _{−j }=(x _{1,−j },…,x _{ N,−j })^{T}. Define \(\boldsymbol {b}^{(m)}=\left (\left (1z_{1j}\right)\left (x_{1j}\mu _{1j}^{(m)}\right), \ldots, \left (1z_{Nj}\right)\left (x_{Nj}\mu _{Nj}^{(m)}\right)\right)\) and \(E^{(m)}=X_{j}^{T}\text {diag}\left (\left (1z_{1j}\right)\mu _{1j}^{(m)}, \ldots, \left (1z_{Nj}\right)\mu _{Nj}^{(m)}\right)X_{j}\). If we ignore additive constants, the quadratic approximation of the objective function at \({\hat{\boldsymbol{\beta}}}_{j}^{(m)}\) yields
for an appropriate σ ^{(m)}. To find an appropriate upper bound, we may set σ ^{(m)} as the maximum of \(\left (1z_{ij}^{(m)}\right)\mu _{ij}^{(m)}\) for i=1,…,N. We can easily show that
is a positive definite matrix. The upper bound can be expressed as
where \(\boldsymbol {w}_{j}^{(m)}=X_{j}\hat {\boldsymbol {\beta }}_{j}^{(m)}+{\sigma ^{(m)}}^{1}{b}^{(m)}\). The majorized problem is written as
Up to a constant depending not on β _{−j } but on \({\hat{\boldsymbol{\beta}}}_{j}^{(m)}\), the function in the minimization problem (4) majorizes \(l_{j}^{c2}\). Hence we achieve the property, guaranteeing the convergence of the algorithm for \(\boldsymbol {\beta }_{j}^{(m)}\) in Mstep.
Results
In this section, we illustrate that our method is effective through a simulation study by comparing the performances of our method, LPGM, and NPGM on simulated data. Then we apply our method to a RNA sequencing data. Also we investigate the gender effect by comparing the estimated networks according to different genders.
Simulation
To simulate data from a Poisson network with excess zeros, we modify the data generation scheme in [4] slightly. Let X∈{0,1,⋯,∞}^{N×p} denote n independent observations from a Poisson network with p nodes. The data generation model is given as
where Y is a N×(p+pC _{2}) matrix with Y _{ ij }∼iidPoisson(λ _{true}) and E is a N×p matrix with E _{ ij }∼iidPoisson(λ _{noise}). The coefficient matrix B encoding the true underlying graph structure denoted by the adjacency matrix A∈{0,1}^{p×p} is defined as
where P is the p×pC _{2} pairwise permutation matrix, ⊙ denotes the elementwise product, and tri(A) is the pC _{2}×1 vectorized upper triangular portion of the adjacency matrix. Each of offdiagonal elements in A is randomly generated from Bernoulli(ρ), where ρ is the sparsity parameter for the network defined as the number of active edges in A divided by the number of all possible edges between the nodes. In order to make X _{ ij }’s to have excess zeros, we multiply each of X _{ ij } by a random variate from Bernoulli(π) for i=1,…,N and j=1,…,p. As an abuse of notation, we denote the final matrix containing zero inflated Poisson counts as X.
The Poisson rates were set as λ _{true}=1.5 and λ _{noise}=0.5. And we have experimented at different levels of N(=50,100,150), p(=10,20,30), π(=0%,10%,20%), and ρ(=.2,.3,.4). At each experimental condition, we generated data according to the above scheme and compared the areas under the curve (AUC) from ZILPGM, LPGM, and NPGM. AUC can be obtained in this way. If we regard active and inactive edges in A as positive and negative examples in a binary classification, then we can compute true positive rate (TPR) as the fraction of edges found by a method that are in the true underlying network structure A. False positive rate (FPR) can obtained analogously. Receiver operating characteristic (ROC) curve and AUC can be obtained from TPR and FPR. To assess the variabilities, we replicated the process of generating data and computing AUC’s 100 times. In Tables 1 and 2, average AUC’s of ZILPGM and LPGM with their standard errors in parentheses and pvalue from the paired sign rank test on AUC’s over 100 replications are reported.
Let us consider the effects of each factor with the other factors held fixed. As ρ increased (or the network became dense), AUC’s of all the compared methods have decreased. Similarly, as the dimension p grew larger, their AUCs became smaller. As the sample size N grows, AUC’s tends to improve. However the tendency is sometimes not so clear. Now consider the effect of excess zeros. When there is no zero inflation (π=0), AUC’s from ZILPGM, LPGM, and NPGM were not significantly different. When we have zero inflations (π=0.1,0.2), ZILPGM seems to significantly outperform LPGM and NPGM. NPGM seems to be outperformed by LPGM. The gaps between AUC’s from ZILPGM and LPGM when π=0.2 was not necessarily larger than that when π=0.1. A potential explanation for this phenomenon follows. As π increases, we have more zero counts in the data and thus the estimation accuracy for the mixing parameter will improve. Meanwhile, the estimation accuracy for the Poisson parameters can degrade because Poisson parameters are learned from nonzero counts. The tradeoff between these two estimation errors may occur at a certain level of π.
Chromosome data
To investigate the validity of the proposed method, we applied it to the RNA sequencing data in the form of a count matrix that contains the number of mapped reads for 60 normal individuals in [25]. We selected 899 genes in the sex chromosomes, i.e., X and Y, first. Each of 899 genes has many zero counts. For a gene with almost all the counts equal to zero, its mixing parameter is estimated as one. To reduce the computation, we have reduced the original data to a data of dimension n=60,p=360 by keeping genes with the number of nonzero counts less than or equal to one.
Figure 1 shows the estimate for the network structure from our method. While 49 genes are clustered together, the other genes remain isolated. Top ranked genes are shown in Table 3 according to their degrees. Note that the degree of a gene is the number of edges being incident upon the gene. We further identified the function of genes with large degree. By GOBP annotation, NDUFA1 and NDUFB11 are involved in mitochondrial electron transport chain (especially complex I), which affects the capacity for the production of ATP through oxidative phosphorylation. GO annotations related to MID1IP1 and PIM2 are protein Cterminus binding and transferase activity, respectively. Proteins with these functions should highly interact with other proteins to control regulation process in cells. Meanwhile, genes with small degrees, SYAP1 and P2RY10, involved in PI3K/Akt signaling and Gprotein coupled receptor (GPCR) activity, respectively [26]. GPCR activate the PI3K/Akt signaling pathway involved in the cellular responses including metabolism, proliferation, apoptosis, and survival [27].
Now let us investigate the effect of gender. In order to compare the networks for different gender groups with 27 males and 33 females, we applied our method to each of gender groups separately. The estimated networks for male and female groups are shown in Figs. 2 and 3. The differentially expressed genes in each group are listed in Table 4. Originally, a differentially expressed gene in a treatment and control groups is a gene with mean expression levels in those groups are significantly different. Although our method is not explicitly related to a hypothesis testing for comparing mean levels, it is implicitly related to a hypothesis testing for conditional independence of counts between genes through a regularized graph learning method on count data. So, in this sense, we call a gene differentially expressed in male and female groups if it appears only one of the the estimated networks for male or female groups. For example, ARMCX1 was selected as a node in the network of the male group and CLIC2 was selected in the network of the female group.
The ARMCX1 gene encodes a member of the ALEX family of proteins and is located on the X chromosome. It was reported that downregulated ARMCX1 transcripts have been found to be significantly reduced prostate cancer and may play a role of tumor suppressor gene [28, 29]. CLIC2, a member of the glutathione Stransferase structural family and a suppressor of cardiac ryanodine receptor (RyR2) Ca2+ channels located in the membrane of the sarcoplasmic reticulum, is controlled by redoxdependent processes and would allow to limit cellular damage in terms of oxidative stress [30]. Above mentioned cellular oxidant detoxification and glutathione metabolic process could inhibit agerelated deterioration, protect the human neuronal cells, and regulate the expression of many genes primarily involved during immune system activities and inflammatory responses [31].
Following GO functional enrichment analysis, genes differentially expressed in the male group included SLC9A7, PLP2, MAGT1, COX7B, STK26, CYBB, MMGT1, BCAP31, and SLC9A6, whereas genes differentially expressed in the female group were RAB33A and UBQLN2. The differentially expressed genes in the male group are involved in the ion transportrelated pathways, whereas the differentially expressed genes in the female group are involved in the regulation of autophagosome assembly.
It has been implicated that ion transport pathways may play a key role in the male reproductive potential, such as capacitation and the acrosome reaction, which are critical steps in sperm physiology preparing for fertilization [32]. On the other hand, it has been investigated on the formation of an autophagosome stimulated by oxidative or metabolic stress taking into account the sex/gender disparities in terms of immunity and inflammation [33–35]. Furthermore, these advantages of women in immunity and inflammation have been well known and these phenotypic differences in immune responses from males result from direct genetic differences [34, 36].
Discussion
In this paper, we propose a penalized version of zero inflated spatial Poisson regression and derive an efficient EM algorithm built on coordinate descent. On simulated data, our method was shown to yield competitive performances in terms of AUC. Particularly, in the presence of excess zeros, our method outperformed LPGM, which is a state of art method in learning graph structures for count data. Note that one may apply the likelihood ratio test for nonnested hypotheses in [37] in order to test for excess zeros on each node. Also we have applied our method to the chromosome data to infer its network structure. Constructing the networks for different genders, we identified the genes differentially expressed in the male and female groups.
There are several issues we have not addressed in this paper. First, one may study the properties our estimators. For Guassian graphical models, asymptotic properties of the estimators are rather well studied in the literature. For example, [38] study asymptotic normality and optimalities in the estimation of Gaussian graphical models. Monod [14] provides sufficient conditions for the consistency of the MLE for ZISP model and discusses some properties such as asymptotic normality and efficiency of the MLE. Because our model is based on a penalization of ZISP model, the results in [14] will provide a starting point for studying properties of the estimators. Particularly, in our case, the properties of the estimators for the incidence matrix rather than the coefficients are of interest. Second, our method can be applied to construct biological networks as well as other networks for count data with excess zeros. Examples include userratings, spatial incidence of a disease or crime, worddocument counts, and others. Third, one may also extend our model to Poisson graphical models with multipleinflations as in [39]. Still another direction is to generalize our model to other distributions such as negative binomial and gamma distributions.
Conclusions
In the present study, expression of ARMCX1 and CLIC2 turned out to be different according to gender. Very little is known about the functional properties of these two genes, this could make ARMCX1 and CLIC2 the possible candidates of medical relevance, such as prostate cancer in male [28, 29] and oxidative stressrelated diseases for female [40]. Therefore, further evidences seem to be necessary for identifying gene expression patterns and validating its diagnostic potential that differentiated patients with relevant diseases from healthy controls in each sex in the populationbased cohorts and, afterwards, it will be translated to clinical practice with its diagnostic impact.
Abbreviations
 AUC:

Area under the curve
 EM:

Expectationmaximization
 FPR:

false positive rate
 GPCR:

Gprotein coupled receptor
 LLA:

Local linear approximation
 LPGM:

Local Poisson graphical model
 MM:

Majorization minimization
 NPGM:

Nonparanormal graphical model
 PGM:

Poisson graphical model
 ROC:

Receiver operating characteristic
 TPR:

True positive rate
 ZILPGM:

Zeroinflated local Poisson graphical model
 ZIP:

Zeroinflated Poisson
 ZISP:

Zeroinflated spatial Poisson
References
 1
Segal E, Wang H, Koller D. Discovering molecular pathways from protein interaction and gene expression data. Bioinformatics. 2003; 19:264–72.
 2
Kotera M, Yamanishi Y, Moriya Y, Kanehisa M, Goto S. Genies: gene network inference engine based on supervised analysis. Nucleic Acids Res. 2012; 40:W162–7. doi:10.1093/nar/gks459.
 3
Gallopin M, Rau A, Florence J. A hierarchical poisson lognormal model for newtork inference from rna sequencing data. PLoS ONE. 2013; 8:431–44.
 4
Allen G, Liu Z. A local poisson graphical model for inferring networks from sequencing data. IEEE Trans NanoBiosci. 2013; 12:1–10.
 5
Besag J. Spatial interaction and the statistical analysis of lattice systems. J R Stat Soc B. 1974; 36:192–236.
 6
Yang E, Ravikumar P, Allen GI, Liu Z. On poisson graphical models In: Welling M, Ghahramani Z, editors. Advances in Neural Information Processing Systems. La Jolla: NIPS Foundation: 2013. p. 1718–26.
 7
žitnik M, Zupan B. Gene network inference by fusing data from diverse distributions. Bioinformatics. 2015; 31:230–9.
 8
She Y, Tang S, Zhang Q. Indirect Gaussian graph learning beyond Gaussianity. 2016. arXiv:1610.02590 [stat.ML].
 9
Liu H, Lafferty J, Wasserman L. The nonparanormal: semiparametric estimation of high dimensional undirected graphs. J Mach Learn Res. 2009; 10:2295–328.
 10
Lambert D. Zeroinflated poisson regression, with application to defects in manufacturing. Technometrics. 1992; 34:1–13.
 11
Buu A, Johnsonb NJ, Li R, Tand X. New variable selection methods for zeroinflated count data with applications to the substance abuse field. Stat Med. 2011; 30:2326–40.
 12
Dobbie MJ, Welsh AH. Modelling correlated zeroinflated count data. Aust NZ J Stat. 2001; 43:431–44.
 13
Mullahy J. Specification and testing of some modified count data models. J Econom. 1986; 33:341–65.
 14
Monod A. A quasilikelihood approach to zeroinflated spatial count data. PhD thesis. Lausanne, Switzerland: École Polyetechnique Fédérale de Lausanne. 2012.
 15
Wang Z, Ma S, Wang CY, Zappitelli M, Devarajan P, Parikh C. Em for regularized zeroinflated regression models with applications to postoperative morbidity after cardiac surgery in children. Stat Med. 2014; 33:5192–208.
 16
Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the em algorithm. J R Stat Soc B. 1977; 39:1–38.
 17
Kim SJ, Dix DJ, Thompson KE, Murrell RN, Schmid JE, Gallagher JE, Rockett JC. Effects of storage, rna extraction, genechip type, and donor sex on gene expression profiling of human whole blood. Clin Chem. 2007; 53:1038–45.
 18
Tian Y, Stamova B, Jickling GC, Liu D, Ander BP, Bushnell C, Zhan X, Davis RR, Verro P, Pevec WC, Hedayati N, Dawson DL, Khoury J, Jauch EC, Pancioli A, Broderick JP, Sharp FR. Effects of gender on gene expression in the blood of ischemic stroke patients. J Cereb Blood Flow Metab. 2012; 32:780–91.
 19
Ronen D, Benvenisty N. Sexdependent gene expression in human pluripotent stem cells. Cell Rep. 2014; 8:923–32.
 20
Siegel C, Turtzo C, McCullough LD. Sex differences in cerebral ischemia: possible molecular mechanisms. J Neurosci Res. 2010; 88:2765–74.
 21
Meinshausen N, Bühlmann P. Stability selection. J R Stat Soc B. 2010; 72:417–73.
 22
Monod A. Random effects modeling and the zeroinflated poisson distribution. Commun Stat Theory Methods. 2014; 43:664–80.
 23
Hunter D, Li R. Variable selection using mm algorithms. Ann Stat. 2005; 33:1617–42.
 24
Lange K. MM Optimization Algorithms. Philadelphia: SIAM; 2016.
 25
Montgomery SB, Sammeth M, GutierrezArcelus M, Lach RP, Ingle C, Nisbett J, Guigo R, Dermitzakis ET. Transcriptome genetics using second generation sequencing in a caucasian population. Nature. 2010; 464:773–7.
 26
GeneCards. http://www.genecards.org/cgibin/carddisp.pl?gene=KMT2D&keywords=KMT2D. Accessed 31 May 2017.
 27
Martelli AM, Evangelisti C, Chiarini F, McCubrey JA. The phosphatidylinositol 3kinase/akt/mtor signaling network as a therapeutic target in acute myelogenous leukemia patients. Oncotarget. 2010; 1:89–103.
 28
Jiang M, Li M, Fu X, Huang Y, Qian H, Sun R, Mao Y, Xie Y, Li Y. Simultaneously detection of genomic and expression alterations in prostate cancer using cdna microarray. Prostate. 2008; 68:1496–509.
 29
Kurochkin IV, Yonemitsu N, Funahashi SI, Nomura H. Alex1, a novel human armadillo repeat protein that is expressed differentially in normal tissues and carcinomas. Biochem Biophy Res Commun. 2001; 280:340–7.
 30
Jalilian C, Gallant EM, Board PG, Dulhunty AF. Redox potential and the response of cardiac ryanodine receptors to clic2, a member of the glutathione stransferase structural family. Antioxid Redox Sign. 2008; 10:1675–86.
 31
Pastore A, Federici G, Bertini E, Piemonte F. Analysis of glutathione: implication in redox and detoxification. Clin Chim Acta. 2003; 333:19–39.
 32
Shukla KK, Mahdi AA, Rajender S. Ion channels in sperm physiology and male fertility and infertility. J Androl. 2012; 33:777–88.
 33
Levine B, Mizushima N, Virgin HW. Autophagy in immunity and inflammation. Nature. 2011; 469:323–35.
 34
Libert C, Dejager L, Pinheiro I. The x chromosome in immune functions: when a chromosome makes the difference. Nat Rev Immunol. 2010; 10:594–604.
 35
Lista P, Straface E, Brunelleschi S, Franconi F, Malorni W. On the role of autophagy in human diseases: a gender perspective. J Cell Mol Med. 2011; 15:1443–57.
 36
Arnold AP. Sex chromosomes and brain gender. Nat Rev Neurosci. 2004; 5:701–8.
 37
Vuong QH. Likelihood ratio tests for model selection and nonnested hyoptheses. J Cereb Blood Flow Metab. 1989; 57:307–33.
 38
Ren Z, Sun T, Zhang C, Zhou HH. Asymptotic normality and optimalities in estimation of large gaussian graphical models. Ann Statist. 2015; 43:991–1026.
 39
Su X, Fan J, Levine RA, Tan X, Tripathi A. Multipleinflated poisson model with l _{1} regularization. Stat Sinica. 2013; 23:1071–90.
 40
Singh H. Two decades with dimorphic chloride intracellular channels (clics). FEBS Lett. 2010; 584:2112–21.
Acknowledgements
Not applicable.
Funding
The research of S. Won was supported by the BioSynergy Research Project (NRF2017M3A9C4065964) of the Ministry of Science, ICT and Future Planning through the National Research Foundation. The research of YJ Kim was supported by the Ministry of Education, Science and Technology of Korea (No. 2012M3A9C4048761). The research of C. Park was supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education (No. 2015R1D1A1A01059984).
Author information
Affiliations
Contributions
Statistical modeling: CP. Algorithm and software development: HC. Data processing and bioinformatics: JG. Biological interpretation: SW, YJK. Simulation and data analysis: SW, SK. Manuscript drafting: CP, HC, YJK. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Changyi Park.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional information
Availability of data and materials
The chromosome data can be downloaded from http://jungle.unige.ch/rnaseq_CEU60/ and all the source codes for the current study are available from the corresponding author.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Choi, H., Gim, J., Won, S. et al. Network analysis for count data with excess zeros. BMC Genet 18, 93 (2017). https://doi.org/10.1186/s128630170561z
Received:
Accepted:
Published:
Keywords
 Count data
 EM algorithm
 Network
 Zero inflation