Bulk gene expression experiments relied on aggregations of thousands of cells to measure the average expression in an organism. Advances in microfluidic and droplet sequencing now permit expression profiling in single cells. This study of cell-to-cell variation reveals that individual cells lack detectable expression of transcripts that appear abundant on a population level, giving rise to zero-inflated expression patterns. To infer gene co-regulatory networks from such data, we propose a multivariate Hurdle model using a finite mixture of singular Gaussian distributions. This permits inference of statistical independences in zero-inflated, semi-continuous data to learn undirected, Markov graphical models. The node-wise conditional log-likelihood of the multivariate Hurdle model is convex and tractable, and allows neighborhood selection. We apply penalized maximum pseudo-likelihood using a group lasso penalty to infer conditional independences. The method is demonstrated in a data set of selected expression in T follicular helper cells, and a high-dimensional profile of mouse dendritic cells. It reveals network structure not present using other methods; or in bulk data sets.