{'fraction_non_alphanumeric': 0.039479942920564455, 'fraction_numerical': 0.013054278315099624, 'mean_word_length': 4.398430813124109, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Nonparametric Multi-view Model for Estimating Cell Type-Specific Gene Regulatory Networks

We present a Bayesian hierarchical multi-view mixture model termed Symphony that simultaneously learns clusters of cells representing cell types and their underlying gene regulatory networks by integrating data from two views: single-cell gene expression data and paired epigenetic data, which is informative of gene-gene interactions. This model improves interpretation of clusters as cell types with similar expression patterns as well as regulatory networks driving expression, by explaining gene-gene covariances with the biological machinery regulating gene expression. We show the theoretical advantages of the multi-view learning approach and present a Variational EM inference procedure. We demonstrate superior performance on both synthetic data and real genomic data with subtypes of peripheral blood cells compared to other methods.


Introduction
Joint analysis of different types of data that are associated with the same underlying phenomenon is more informative than analysis of individual data types, and increases signal to noise ratio (Xu et al., 2013;Wang et al., 2013;Rey & Roth, 2012). This approach, known as multi-view learning or learning with multiple distinct features, has been successfully used in various settings (Li et al., 2002;Jones & Viola, 2003;Hardoon et al., 2004;Pan et al., 2007).
Here, we apply such a multi-view learning approach to address an important biological problem by integrating two views of gene regulation. Our goal is to infer cell clusters (characterizing cell types) as well as their underlying gene regulatory networks (GRNs), which are directed weighted networks between genes depicting the extent to which a regulatory gene influences the expression of each of its downstream target genes. Understanding differences between regulatory mechanisms across different cell types provides valuable insight in normal development of cell types (Davidson, 2010), and mechanisms disrupted in cancer cells (Pe'er & Hacohen, 2011;Kreeger & Lauffenburger, 2009).
Recent advances in single-cell genomic technologies (Hashimshony et al., 2012;Jaitin et al., 2014;Shalek et al., 2013) which measure gene expression at the resolution of individual cells, present remarkable opportunities to characterize different cell types by clustering cells based on heterogeneity of gene expression (as observed features) Macosko et al., 2015). Learning GRNs from gene expression data alone, however, leads to detection of spurious network links based on correlated genes, while integrative learning from multiple data sources has been shown to improve overall joint inference (Zhu et al., 2008;Hecker et al., 2009;Azizi et al., 2014). Therefore, we aim to identify GRNs driving heterogeneous cell types through integrating single-cell expression data with other genomic data types. In particular, epigenetic technologies such as ATAC-seq (Buenrostro et al., 2015a), scan the genome for accessible DNA regions, identifying potential interaction between a gene and regulator proteins translated from other genes. In other words, epigenetic data contains information about direct regulatory links between genes and incorporation of epigenetic data is a promising direction for improved inference of GRNs (Guo et al., 2017;Rotem et al., 2015).
We present a novel integrative model, which we refer to as Symphony, as a Dirichlet process mixture model that jointly learns clusters of cells and GRNs specific to each cluster. Symphony is an extension of the BISCUIT model (Prabhakaran et al., 2016;Azizi et al., 2018) which clusters cells while simultaneously distinguishing biological heterogeneity from technical noise in single-cell gene expression data. This is done through incorporating cell-specific parameters scaling the cluster means and covariances for a multivariate Gaussian mixture model.
We extend the BISCUIT model and replace the hyperparameters with a generative process exclusively driven by the paired epigenetic data, which captures the biological mech-arXiv:1902.08138v1 [stat.ML] 21 Feb 2019 Figure 1. The generative process in Symphony: We aim to infer gene regulatory networks (GRNs) denoted as R k specific to each cluster of cells from integration of two views: Epigenetic data (indicative of network edges) from the bulk of cells (C) and expression of genes (network nodes) at resolution of single cells (X). GRNs are directed weighted networks with edge weights denoting regulatory impact of one gene on another (e.g. activation or repression); the impact of regulation is reflected in covariance between genes (Σ k ).
anism responsible for observed gene covariances per cell type. Briefly, the epigenetic profiles which denote accessible DNA in the bulk samples are deconvolved into cell-type specific accessible regions (Figure 1). Within these regions, the binding of regulatory proteins translated from genes impacts the expression of nearby genes, such that accessible regions may be mapped to gene-gene interactions. This mapping is based on prior knowledge of recurring DNA sequences (known as motifs) associated with these regulatory proteins which occur in regions of accessible DNA. Most importantly, the covariance in observed gene expression is related to a graph power of the regulatory network, capturing the propagated impact of regulation in the network (indirect regulation)( Figure 2). This multi-view framework can also be applied to other settings, such as text characterization. For example, to learn the context of queries (vector of words), the bag-ofwords simplification may not be sufficient (Biemann, 2005). However, the order of words, which can be represented as a latent directed network can imply the context, and incorporating observations from this network such as the frequency of one word following another (as a second view) can enhance extraction of context and clustering of queries (observed as the first view) (Landauer et al., 1997;Recchia & Jones, 2009).
Related technologies and methods. The problem of inferring GRNs specific to cell types involves identifying differences in gene-gene interactions across cell types. However, in most cases the cell types are not well-characterized, hence gene markers are not known to enable sorting of cell types prior to measuring epigenetic data (gene-gene interactions). Therefore, epigenetic data measured on the bulk of cells represents a mixture of cell type-specific epigenetic profiles. One solution is measuring epigenetic profiles at the resolution of single cells. These technologies have only recently emerged (Buenrostro et al., 2015b); we therefore constructed a model to allow integration of bulk ATAC-seq data. Symphony can be easily adapted for inferring GRNs from single-cell ATAC-seq data as well.
Other works have attempted to apply computational deconvolution algorithms intended for bulk expression data, such as those using source separation techniques (Houseman et al., 2016), NMF-based methods (Repsilber et al., 2010) or Bayesian models (Erkkilä et al., 2010), to instead infer cell type-specific epigenetic profiles. Recent methods such as SCENIC (Aibar et al., 2017) infer GRNs from single-cell expression data alone and do not incorporate epigenetic or other types of data.
Contributions. In this paper, we show that a multi-view learning framework would improve the deconvolution of epigenetic data. Furthermore, using an integrative model, we improve the clustering of cells, and hence characterization of cell types. Most importantly, our model presents the advantage of inferring cell type-specific GRNs that give insight into heteroegeneity of underlying mechanisms across cell types. We present a Variational EM inference procedure and show that the integration guarantees model identifiability, while learning from the epigenetic view alone does not. While other works have attempted to integrate bulk multiomics data (Lake et al., 2017;Brown et al., 2013;Ritchie Figure 2. Symphony captures direct and indirect regulation. The impact of regulation (R k ) is propagated through the network up to path length of two and is reflected in covariance between indirectly connected genes (Σ k ).
Figure 3. Validation of model assumption. Weighted sum of actual peak heightsp p p k s measured from sorted clusters of CD34+ hematopoeitic cells using ATAC-seq, with weights proportional to proportions of cell types, compared to measured peak height from ATAC-seq on the bulk of cells c c c; heatmap shows density with yellow (low) to blue (high). et al., 2015), there are no methods to our knowledge that infer heterogeneous GRNs through integrating epigenetic and single-cell resolution gene expression data.

Model
The observed data is considered as two views from the biological system ( Figure 1): View 1. View 1. View 1. Single-cell gene expression data from scRNA-seq technologies (Klein et al., 2015;Macosko et al., 2015) denoted as X d×n = [x x x 1 , · · · , x x x j , · · · , x x x n ] where each observation x x x j ∈ R d for cell j ∈ {1, · · · , n} corresponds to d genes (as features). Each entry x ij for i = [1, · · · , d] contains the expression of gene i in cell j (more precisely, the log of counts of mRNA molecules per gene i from cell j plus a pseudo-count).
View 2. View 2. View 2. Epigenetic data, for example measured with ATAC- seq technology (Buenrostro et al., 2015a), denoted as C l×r = [c c c 1 , · · · , c c c t , · · · , c c c r ] where each observation c c c t ∈ R l for t ∈ {1, · · · , r} corresponds to l genomic regions (as features). Specifically, c c c t is an experimental replicate measuring accessibility of genomic regions m = [1, · · · , l].
Prior knowledge. The genomic regions in C can be mapped to genes in X with a pre-defined mapping function g(i, i ) = m that relates each genomic region m ∈ {1, · · · , l} to a gene-gene interaction i → i for i, i ∈ {1, · · · , d}. We also define M d×d based on prior knowledge containing binary values M i,i = 1 if the motif sequence for gene i exists in the genomic region m in the vicinity of gene i, meaning a potential interaction can exist from gene i to gene i.

Epigenetic Model (View 2)
The epigenetic data is informative of network structure, i.e. existence of edges between genes (features) and gene expression data contains information on network nodes. We aim to infer this directed weighted network (GRN) for each cluster k ∈ {1, · · · , K} of cells, denoted by the asymmetric matrix R d×d k in which entry R k i,i = 0 if i → i, meaning gene i directly regulates gene i. R k i,i is the regulatory function of gene i on gene i in cluster k such that R k i,i > 0 or R k i,i < 0 represent activation or repression of expression respectively, with |R k i,i | being the strength of regulation.
We do not aim to distinguish all layers of the regulatory process (such as protein phosphorylation) and rather interpret GRNs as an approximation for the overall impact of TFs on target genes at the transcriptional level.
We model regulation of gene expression as follows: Genome accessibility in cluster k is represented with latent variable p p p k = [p 1 k , · · · , p l k ] T ∈ R + l containing l genomic regions (features). This represents log of peak heights plus 1 (to ensure positive domain) in all genomic regions, for each cell type k. We set a truncated multivariate Normal prior to capture the structure between genomic regions encompassing co-regulated genes (i.e. genes sharing regulators) with mean η η η and covariance Λ: p p p k ∼ truncN (η η η, Λ, 0 0 0, +∞). In this paper, we assume a setting where we do not observe p p p k s, and only observe epigenetic data from the bulk of cells which can be represented as a weighted sum of clusterspecific epigenetic profiles where π k s are weights. Thus, our epigenetic model is: We validated the above assumption of weighted sum using ATAC-seq data from hematopoeitic progenitor cells from Corces et al. (2016) by computing the weighted sum of measurements from sorted cell types (actual p p p k s denoted withp p p k ) to measurements from the bulk of cells as c c c ( Figure  3).
We have a K-order Dirichlet prior over π k : π k |ϕ, K ∼ Dir (π k | ϕ K , · · · , ϕ K ), where ϕ −1 ∼ Gamma(1, 1). Then, in cell type k, if a genomic region m in the vicinity of gene i is accessible with log peak height p m k , and the motif sequence associated with one or more Transcription Factor (TF) proteins translated from genes i , i , ... exists in the region (M i,i = M i,i = 1), then the TF(s) can bind to the region and hence regulate the expression of gene i. Furthermore, we assume the peak height (strength of TF binding to genome) is informative of i → i edge weight |R k | (strength of regulation). Thus, we model R k as follows: The function g(·) maps gene pair i, i to genomic region l. S denotes a sign indicator variable representing repression or activation function. We set S according to the sign of the empirical covariance: S i,i = sign(Σ i,i ).

Single-cell Gene Expression Model (View 1)
We use the above epigenetic model to drive gene expression data based on the following key ideas: First, if a direct regulatory link exists from a TF associated with gene i to a target gene i (i → i), we assume there is strong covariance between their expressions. Second, due to van der Corput's inequality (Montgomery, 2001), covariances can partly reflect the propagated impact of indirect regulation in cases where genes are not directly connected but exist on the same path in the network ( Figure 2). For example if i → i and i → i, we might also observe covariance between i, i even though they are not directly connected in the network (e.g. Here, we consider indirect effects with path length up to two using the square of the indirected network (R k + R T k ) 2 (Walker, 1992) such that: 2 is positive semi-definite according to Lemma 3 in the following section making the above modeling assumption feasible. Additionally, this model can capture combinatorial regulation in the inferred covariances. In particular, a gene pair i, i will always have the same directionality of regulation (i.e. activation or repression relationship), but Σ i,i k can be positive in one cluster and negative in another cluster depending on the relative regulatory strength of activators and regulators in its path. An example of this variability in sign is shown in Supplmentary Figure 12.
Gene expression data for each cell j denoted as x x x j is then modeled similar to the multivariate Gaussian mixture model in BISCUIT: where z j denotes assignment of cell j to cluster k ∈ {1, · · · , K}. With integrating two data modalities (gene expression and epigenetic data), we improve inference of clusters (cell types). The scaling parameters α j , β j specific to cell j are used to normalize the data x x x j in downstream analysis by transforming to y y y j ∼ N (µ µ µ k , Σ k ) according to the cluster it is assigned to z j = k, similar to Prabhakaran et al. (2016). The plate model for Symphony is summarized in Figure 4.

Theory
We show theoretical advantages of integration of the two data types using Symphony as follows: We define f (x x x) := N (αµ µ µ k , βΣ k ) ∈ R d as the multivariate Gaussian density of x x x and f (c c c) := N ( k π k p p p k , ζI) ∈ R l as the multivariate Gaussian density of c c c. First, we emphasize that the epigenetic model alone f (c c c|p p p k , π k , ζ) is not identifiable and therefore precise inference of deconvolved epigenetic profiles (p p p k s) is not possible: This motivated us to build an integrative model. Identifiability of the single-cell expression model f (x x x|µ µ µ k , Σ k , α α α, β β β, z z z) has been shown under certain enforced constraints on both α, β α, β α, β: While the above conditions guarantee identifiability, they are not inferred from data or biologically motivated and hence interpretation of parameters may not provide the best characterization of cell types. Here, we show that in the integrative model, the constraints for β j s are no longer required and identifiability of the expression model f (x x x|µ µ µ k , Σ k , α α α, β β β, z z z, p p p k , R k ) is guaranteed through the extension of the model that captures regulation, from which we observe additional epigenetic data c c c. Hence the full model f (X, C|{µ µ µ k , Σ k , α α α, β β β, z z z, p p p k , R k , π k , ζ}) is identifiable. We will use the following lemma to show (R k + R T k ) 2 is positive semi-definite.
Lemma 3 Square of a symmetric matrix H gives a symmetric positive semi-definite matrix L (Proof in Supplementary section B). Since single-cell expression data X usually contains expression of genes that do not have observations in their mapped genomic regions in c c c such that d > l, we first define a reduced version of the model where all genes in x x x do have mapped genomic regions: being subsets of µ µ µ and Σ. We next show the identifiability of the reduced model. Then, we use this result to extend the identifiability to the full model given β β β which is the parameter scaling Σ k . Finally, we show the identifiability of the full model.

Inference
We applied the co-ordinate ascent mean field variational inference (CAVI) (Blei et al., 2017;Ghahramani & Beal, 2001) which assigns independent factors to the latent variables. The blueprint for corresponding CAVI updates are below and full derivations are presented in Supplementary section C.
where c is the integration constant (details are included in Supplementary section C). Since the E-step takes ≈ O(d 3 ) due to three matrix inversions, we substitute this with z M AP (x x x) = arg max z p(x x x|z z z)p(z z z|π π π) which takes ≈ O(d 2 ) showing global performance of Symphony in deconvolving epigenetic data; shown with projection of inferred peak heights on principal components of ground truth peak heights (from ATAC-seq on sorted cell types in Corces et al. (2016)) (left) compared to deconvolution using Deconf (right) which fails to deconvolve the majority of peaks shown with overlapping squares (d) Scatterplot of inferred peak heights for all clusters vs ground truth peak height using Symphony (left) compared to Deconf (right); peaks are colored by density; r values show Pearson correlation; peaks constrained by expression data and bulk epigenetic data are triangular points with black outline; NK cells were not included in this plot due to the small cell proportion (< 5%), making deconvolution impossible.
when Σ −1 k s and Σ are apriori Cholesky decomposed. Given the complexity of the model (refer to Supplementary section C), we implemented Symphony using the probabilistic programming language Stan (Carpenter et al., 2016). Furthermore, for a scalable implementation applicable to real genomic data containing thousands of cells, we used the probabilistic programming language, Edward (Tran et al., 2016; 2017) for variational EM (details and approximations presented in Supplementary section E).

Synthetic Data
We first evaluated the performance in deconvolving epigenetic data, clustering cells, and inferring GRNs using data simulated from the Symphony model. We simulated data for n = 100 cells in K = 3 clusters with d ranging from 5 to 20 genes and l = 50 using the Symphony model.
Inference of p p p k , R k . Figure 5(a) shows scatterplots of deconvolved peak heights (p p p k ) compared to actual data. We compared the performance of Symphony to two other deconvolution methods: deconf (Repsilber et al., 2010) which uses NMF, andDsection (Erkkilä et al., 2010) which is based on a Bayesian model. While Dsection captures only the largest cluster, deconf underestimates the cluster-specific peak heights. The behavior of Dsection was reproducible across simulations, and is likely due to the lack of identifiability in the model for epigenetic data alone, as discussed above. Figure 6 summarizes the error in estimating p p p k s across 10 synthetic datasets with the same size as above. This shows the value of incorporating expression data (view 1) in deconvolution of epigenetic data (view 2). Figure 6 also shows a heatmap of inferred R k in one of the synthetic experiments as an example, compared to the actual R k , confirming the ability of Symphony in inferring GRNs.
Clustering performance. We then show the performance in clustering with integrating both views as compared to only using gene expression data (view 1) by computing Fscores across 10 experiments with the same size as above. We compared the performance to BISCUIT (Prabhakaran et al., 2016), as well as other methods commonly used for clustering cells in single-cell gene expression data including DBscan , Phenograph (Levine et al., 2015), Spectral clustering (Ng et al., 2002) and k-means (with K = 3) ( Figure 6). These results show improvement over BISCUIT due the epigenetic extension of the model and significant improvement over other methods. DBscan was unable to cluster the majority of cells, likely due to the small dimensionality of the feature space used in simulations. This shows the value of incorporating epigenetic data (view 2) in improving clustering performance, as compared to using expression data (view 1) alone. This has further value in biological interpretation of clusters as cell types that have both similar expression and similar underlying mechanisms driving expression.

Genomic Data
We also evaluated the performance of Symphony on real genomic data. We used previously published single-cell expression data for peripheral blood mononuclear cells (PBMCs) from Zheng et al. (2017) combined with ATACseq data for PBMCs from Corces et al. (2016). For singlecell expression data, we chose a subset of PBMCs from (Zheng et al., 2017) as X containing n = 6825 cells which express known gene markers for either monocytes, B cells or T cells and NK cells. We chose to focus on d = 28 transcription factors which showed high standard deviation in expression data and are known to be lineage-defining factors. For epigenetic data matching the above cell types, we generated r = 1 mixture of epigenetic measurements with l = 1053 peaks from real ATAC-seq data collected from the above sorted cell types in Corces et al. (2016), with weights π k proportional to frequency of cell types in blood, and used this as observed epigenetic data C. We fixed the clustering in this experiment using Phenograph-derived assignments which we mapped onto nearest cell types to match with epigenetic data (Levine et al., 2015).
We determined non-zero entries in M from ATAC-seq using the FIMO algorithm (Grant et al., 2011), which scans the sequence under the ATAC-seq peak for the occurrence of a motif. We associated a peak with the target gene closest in genomic distance to the peak in these experiments. This assignment is independent of the model structure and can be manually defined by the user.
In the following tests, we used a scalable implementation with Edward (Tran et al., 2016) detailed in Supplementary section E. Figure 14 shows the performance of this imple- Cell type characterization. In this test, we pre-imputed and normalized data, and fix BISCUIT-derived normalization parameters (α, β). Specifically, we normalized and imputed x x x j for each cell j by transforming it to y y y j with y y y j = Ax x x j + b, and setting A = I/β j and b = (I − α j A)µ µ µ k with BISCUIT -inferred parameters. This transformation corrects for cell-specific technical effects captured by α j .β j , as y y y j ∼ N (µ µ µ k , Σ k ) while x x x j ∼ N (α j µ µ µ k , β j Σ k ) (Figure 4) (Prabhakaran et al., 2016). Figure 7 (a) shows t-SNE projections (Maaten & Hinton, 2008) of y y y j s after normalizing expression data based on inferred parameters µ µ µ k , Σ k , α α α, β β β where cells are colored by cluster assignments z j s for K = 5 distinct cell populations. Figure 7 (b) shows normalized expression of known marker genes used to characterize the clusters as monocyte, B, T and NK cell types.
Deconvolving epigenetic data. Figure 7 (c,d) show inferred peak heights p p p k for all clusters using Symphony compared to ground truth cell type specific peak heights (measured with ATAC-seq from sorted cell types). Supplementary Figure 13 shows an example genomic region with differential peaks for three of these cell types showing distinct epigenetic profiles. Projection of peak heights to principal components of ground truth peaks (excluding peaks un-constrained by expression data and peaks which show 0 accessibility in some cell types, to show performance at deciphering magnitudes) shows superior performance in deconvolving all subsets of cells except for the smallest population (NK cells, < 5% of cells). The small error between estimated p p p k s and ground truth cell type-specific peak heights confirms that our model is a good fit for the biological mechanism of regulation. We also evaluated the deconvolution of epigenetic data using deconf (Figure 7 (c)) Figure 11. Identification of cell type-specific regulators. Violinplots of inferred cell type-specific regulatory function R (green) for TFs that show variability across cell types, compared to inferred covariance Σ (pink) between example TFs (regulators) and all of their target genes; positive and negative signs in R denote activation and repression of expression, respectively; values have been z-score normalized to display both variables on the same y-axis.
which shows under-estimation or inaccurate estimation of peak heights. We could not test Dsection as the number of clusters exceeds the number of replicates.
Inferring GRNs. The main advantage of Symphony is the inference of GRNs. Figure 8 shows Symphony can successfully learn cell type-specific covariances comparable to empirical covariances. The gene-gene covariance matrix is then explained by direct regulation as well as the propagated impact of regulation through (R + R T ) 2 , which we previously validated is capable of capturing epigenetic information through the accuracy of its prior p. Furthermore, Figure 9 reveals that inferred regulatory functions are explained by either TF binding strength (peak height) or TF-gene covariance or both. Figure 10 shows the inferred GRNs between TFs with strongest inferred links (|R k | > 4.5) in each cluster. The differences in the structure of the networks suggests different mechanisms driving cell type-specific expression.
We observe numerous regulatory interactions that are variable across clusters. Figure 11 shows examples of TFgene interactions (R) that are also supported by known literature. It can be seen that regulatory functions are partially supported by gene-gene covariances (Σ). We observe CEBP A, CEBP B differentially regulating target genes in monocytes (cluster 5). A recent study (Jaitin et al., 2016) has shown that knock-outs of CEBP B block monocyte differentiation.
We observe regulatory edges in the GRN for T cells between Transcription Factors (TFs) GAT A3, RU N X2 and their target genes, with minimal interaction in B cells and NK cells (Figure 11), and indeed these TFs are known to be associated with activation of cytotoxic T cells (Pearce et al., 2003) and CD8 T cell development (Woolf et al., 2003). We also observe cases such as F OS, F OSB with different regulatory functions despite belonging to the same TF family, showing an example of how expression-derived information can further distinguish genetic interactions which cannot be immediately deciphered from epigenetic data.

Conclusion
We present a hierarchical Bayesian mixture model named Symphony that infers clusters of cells representative of cell types and gene regulatory networks (GRNs) specific to cell types. This is done by modeling the regulatory mechanism driving gene expression in each cell type, and assuming two observations as two views from the system: epigenetic measurements, which are informative of network edges and single-cell gene expression data, informative of network node activity. To the best of our knowledge, this is the first computational method that integrates single-cell expression data with epigenetic data. We provide theoretical justifications for the model and an EM-VI procedure. Symphony shows great performance in clustering cells, deconvolving epigenetic profiles and inferring GRNs in both synthetic and real data from peripheral blood cells and shows superiority to other methods that only address one of these problems. Further, Symphony was able to deconvolve epigenetic data when only one replicate was available through integration of expression data, a potentially common task which would be challenging for any source separation technique. Future iterations of the experiments will allow a Symphony -derived clustering of PBMCs to improve the mapping of cells to cell types, particularly for more similar cell types such as CD4+ and CD8+ T cells and when few genes are considered. Applied to the growing single-cell datasets, Symphony can reveal cell type-specific regulation in normal cells as well as disrupted regulation in cancerous cells.  Figure 12. Examples of synthetic Σ k simulated from the same S matrix for clusters k = 1, 2, 3, showing variability in sign of Σ k that can capture the impact of combinatorial regulations.

B. Extended Theory
Lemma 1 The epigenetic model f (c c c|p p p k , π k ) is non-identifiable Proof sketch. Due to having Kl + K unknowns in the mean parameters while having an l dimensional Normal distribution, where K is the maximum number of allowed clusters, we have an under-determined problem. Thus, we can provide multiple parameter sets π k and p p p k leading to the same Normal distribution for c c c.
Lemma 3 Square of a symmetric matrix H gives a symmetric positive semi-definite matrix L. Proof. We show that there exists a symmetric positive semi-definite matrix L iff there exists a symmetrix matrix H that satisfies H 2 = L. H, L ∈ R d×d . Orthogonal diagonalization of L gives QLQ −1 = D 2 where Q is the orthogonal matrix and D is the square root diagonal matrix diag(l 1 2 1 , · · · , l 1 2 d ). If there exists a matrix H where H = Q * diag(l 1 2 1 , · · · , l 1 2 d ) * Q −1 = QDQ −1 , then QHQ −1 = D. We now write: showing H 2 = L. Next assume H is symmetric and therefore all its eigenvalues are real. For some eigenvalue h of H, l is an eigenvalue of L iff l = h 2 implying all eigenvalues of L ∈ R * where R * = {0} ∪ R + . Further, L is symmetric given H is symmetric. We now have L as symmetric and has non-negative eigenvalues proving that L is positive semi-definite.
Proof sketch. Using Lemma 2 and focusing on the marginal distribution f (X 0 |β β β, µ µ µ 0 k , Σ 0 k , α α α, z z z, π k ) we know that β j Σ 0 k and π k are identified for all j and k. Hence, β 1/2 j Σπ k (Σ 0 k ) 1/2 is identified. Using Lemma 3, (R k + R T k ) 2 is non-negative semi definite which allows us to define (Σ 0 k ) 1/2 from the Wishart distribution. Therefore, β 1/2 j Σπ k (R k + R T k ) is identifiable. We also know that Σπ k p p p k is identified through the marginal distribution for C. Thus, putting the above together, it is not possible to have multiple values for β j since that would require the sum Σπ k (R k + R T k ) to be different for two sets of parameters while each element of this sum can be written as elements of Σπ k p p p k which is identified.
It remains to show identifiability for parameters p p p k , R k , ζ from integrated model on C. ζ is identified due to Normality of C. Identifiability of p p p k will lead to identifiability of R k . We focus on proving identifiability of p p p k .
Using the identifiability of Σ k we have d(d − 1)/2 equations based on generation of Σ k based on p p p k . Considering that p p p k has l unknowns and the relationships are at most polynomials of degree 2, as long as d(d − 1)/2 > l 2 we have an overdetermined system of equations to identify p p p k from Σ k .

C. Variational Inference update equation derivations C.1. Joint distribution for X and C
We use the graphical model in Figure 4 to write down the variational inference equations for Symphony. Note that this is constructed based on conditionally-conjugate priors for µ µ µ and Σ and on space discretised by z z z = {z 1 , · · · , z n } where Figure 14. Performance in deconvolving epigenetic data with Edward implementation. Estimated peak heights p p p k using Symphony for n = 4000 simulated cells in K = 3 clusters with d = 100 genes and l = 550 versus true peak heights, compared to two other deconvolution methods: deconf (Repsilber et al., 2010) and Dsection (Erkkilä et al., 2010); each dot represents a genomic region; heatmap shows density. z j = {j}, j = [1, · · · , k]. The joint is written based on the Markov blanket for each parameter.

C.2. Variational distributions
We now write the factorized distribution q which will approximate the joint distribution in Equation 9. We follow Chapter 10 of Bishop (Bishop, 2006) and the Matrix cookbook Section 8.2 (Petersen et al.).
where ψ is the digamma function and ψ(t) = d d0 ln(Γ(t)) = Γ (t) Γ(t) . Runtime complexity∼ O(1) Therefore, We require that q * (z z z|X) is normalised and that for every observation j, there is only one non-zero z jk ∀ k ∈ {1, · · · , K}. Therefore it is sufficient to normalise each ∆ jk as The expectation for the discrete distribution q * (z jk ) gives the responsibilities r jk for point x j with the current k th cluster's parameters : with the overall runtime complexity for calculating r jk is O(d 3 ).
MAP estimate for z Since the E-step has a runtime ∼ O(d 3 ) due to three matrix inversions, we substitute this as z M AP (x x x) = arg max z p(x x x|z z z)p(z z z|π π π) which has a runtime of ∼ O(d 2 ) when Σ −1 k s and Σ are apriori Cholesky decomposed.
(by taking terms in µ k and µ ) 3. Let us denote R * k = (R k + R T k ) 2 . By construction of R k , R * k is a Gram matrix and a valid scale matrix for Wishart parametrisation.
(by replacing α * j with ln α j and making this a logNormal pdf )

Derivation is similar to that of Equation 38
.
(by replacing β * j with ln β j and making this a logNormal pdf ) 6.
(taking terms only in R k ) 9.

E.3. Approximations
Inference of gene-gene GRNs and covariance matrices are the main goals of Symphony, yet accurate inference of such large matrices involves a number of computational challenges. In particular, constraints on covariance matrices of a multivariate normal distribution are difficult to enforce in the optimization setting of variational inference. For instance, large sparse matrices may very easily become non-singular during optimization, leading to un-defined loss functions.
We use several techniques to solve this problem. For one, we define the Wishart distribution in Edward using the Bartlett Decomposition, rather than the built-in Wishart function of tensorflow, which allows us to more easily define variational parameters. Specifically, we replace the sampling of covariance matrices Σ k ∼ W ishart with a generative model constructing only univariate chi-squared distributions c and normal distributions n, which can be shown to produce a valid sample from the Wishart distribution (Kshirsagar, 1959). In this setting, we define variational distributions corresponding to the dummy variables n and c, as opposed to defining a matrix variate distribution which, during the course of optimization, must fit all the constraints of valid covariance matrices. Initialization of these parameters to large values additionally avoids problems with singularity in most cases.
In addition to the use of the Bartlett Decomposition, the Edward version of Symphony replaces the standard Wishart with a scaled Wishart for added flexibility of the model in the variational inference case. The scaled Wishart necessitates addition of a latent parameter δ, such that Addition of the normal distribution above to the generative process infuses flexibility to the Wishart, whose variance is usually defined by a single degrees of freedom parameter. In addition, we allow separate inference of the diagonal and off-diagonal of the covariance matrices. This is a desirable property for Symphony, in that the model of gene regulation does not necessarily capture the diagonal of the covariance matrices representing variances of gene expression. Likewise, we solve additional issues caused by matrix inversion by simply replacing the prior on Σ k with a Wishart instead of Inverse-Wishart distribution. We note that, while this choice is not conjugate, this is valid as both distributions satisfy the requirements for priors on the covariance matrix.
We require a variational EM procedure with Edward, such that the cluster assignments z are updated every several iterations with a maximization step. In particular, we choose z i for each cell based on the maximum likelihood of cluster assignment in the Gaussian mixture. This prevents the need for discrete optimization over categorical variational parameters. The performance of the variational EM algorithm is maximized when a good initial value for clustering is chosen. In this work, we initialized clusters based on cell-cell kNN graphs with Phenograph (Levine et al., 2015).
Finally, we made some other small distributional changes which seemed to produce better results in our experiments with this particular implementation. We replace the multivariate prior on p with a univariate prior centered at a constant. In addition, we treated binary M as a latent variable with a very tight variance. This was to add additional flexibility to the model, and to also assist with singularity issues by providing a dense matrix mean to R. We note that all fitted values of M were similar within a small to their previously fixed values of either 0 or 1.