Network analysis for count data with excess zeros

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 expectation-maximization (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][2][3][4]. For more details on Markov network inference, see those and the references therein.
*Correspondence: cpark463@gmail.com 6 Department of Statistics, University of Seoul, 02504, Seoul, Korea Full list of author information is available at the end of the article The main focus of this study is to infer the network structure for a count data. The auto-Poisson 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 auto-Poisson 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 zero-inflated 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 well-studied 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 one-step 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 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 X = (x ij ) N×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 μ j = exp β j + k =j β jk x k , β 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 pairwise-only dependency models, the situation is basically similar. In order to have a valid joint distribution, the coefficient for the interaction term β jk should be non-positive. 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 log-likelihood of (1) in the general weighted LASSO form: where μ ij = exp β j + k =j β jk x ik , λ ≥ 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 equal-spaced grid points in log scale on [λ max , λ min ], where λ max = max j∈{1,··· ,p} max k =j 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 log-likelihood 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:Ê = {(j, k) : β jk = 0 orβ kj = 0, j = k}. We devise an EM algorithm as in [15] to minimize (3).

Computational algorithm
Let O j = {i : x ij = 0} and P j = {i : x ij = 0}. The negative log-likelihood function in (2) is the sum of However, it is difficult to maximize this likelihood directly because the score function of − i∈O j log π j / 1 − π j + e −μ ij 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 The log-likelihood 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 (m) −j and the probability of zero state at mth step as Our EM algorithm alternates the following steps until convergence.
• E-step: Estimate z ij by its conditional mean z (m) ij given data and parameters from the previous step.
Here we set the initial values for our EM iteration as π (0) j = the number of zeros of jth variable/n for j = 1, · · · , p and β For each variable, we use the Majorize-Minimization (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 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 c2 If we ignore additive constants, the quadratic approximation of the objective function atβ (m) −j yields for an appropriate σ (m) . To find an appropriate upper bound, we may set σ (m) as the maximum of We can easily show that is a positive definite matrix. The upper bound can be expressed as The majorized problem is written as Up to a constant depending not on β −j but onβ (m) −j , the function in the minimization problem (4) majorizes l c2 j . Hence we achieve the property, guaranteeing the convergence of the algorithm for β (m) −j in M-step.

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 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 element-wise product, and tri(A) is the pC 2 × 1 vectorized upper triangular portion of the adjacency matrix. Each of off-diagonal 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 in-active 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 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 non-zero 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 GO-BP 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 C-terminus 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 G-protein 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 S-transferase structural family and a suppressor of cardiac ryanodine receptor (RyR2) Ca2+ channels located in the membrane of the sarcoplasmic reticulum, is controlled by redox-dependent 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 age-related 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 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][34][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 non-nested 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 user-ratings, spatial incidence of a disease or crime, word-document counts, and others. Third, one may also extend our model to Poisson graphical models with multiple-inflations 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 stress-related 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 population-based cohorts and, afterwards, it will be translated to clinical practice with its diagnostic impact.