Gene Regulatory Network Inference from Single-Cell Data Using Multivariate Information Measures

Abstract

While single-cell gene expression experiments present new challenges for data processing, the cell-to-cell variability observed also reveals statistical relationships that can be used by information theory. Here, we use multivariate information theory to explore the statistical dependencies between triplets of genes in single-cell gene expression datasets. We develop PIDC, a fast, efficient algorithm that uses partial information decomposition (PID) to identify regulatory relationships between genes. We thoroughly evaluate the performance of our algorithm and demonstrate that the higher-order information captured by PIDC allows it to outperform pairwise mutual information-based algorithms when recovering true relationships present in simulated data. We also infer gene regulatory networks from three experimental single-cell datasets and illustrate how network context, choices made during analysis, and sources of variability affect network inference. PIDC tutorials and open-source software for estimating PID are available. PIDC should facilitate the identification of putative functional relationships and mechanistic hypotheses from single-cell transcriptomic data.

Keywords: gene regulation; mutual information; network reconstruction; single-cell PCR; single-cell RNA-seq.

Figure 1

Network Inference from Single-Cell Data Single-cell transcriptomic data quantify mRNA species present inside individual cells. By considering pairs (or triplets, quadruplets, etc.) of mRNA species, we can test for statistical relationships among them. These dependencies may reflect coordinated gene expression of these pairs (or groups) of genes, resulting from gene regulatory interactions or co-regulation. Once such sets of genes that jointly change in expression are known, other statistical, bioinformatic, or text-mining analyses can be used to identify likely transcriptional regulators for these sets of genes. By iterating such in silico analyses with further, targeted experimental studies, we can, in principle, build up a representation of the gene regulatory network.

Figure 2

PID Profiles for Three-Gene Networks (A) Mean PID values for three-gene networks with different topologies. PID values are calculated using data simulated from three-gene networks with the topologies illustrated above each plot; the models used for simulation assumed mass action (green) or thermodynamic (orange) kinetics. For each three-gene network, 12 PID values were calculated from the simulated data — there are four PID contributions with each gene treated as the target gene in turn, consisting of a redundant, synergistic, and two unique contributions (Equation 5 in Box 1). Each line graph shows the mean PID values calculated from simulations using different initial conditions (error bars indicate 1 SD), with the horizontal axis labels indicating the PID contribution, e.g., the first four values show the PID values with gene X as the target, consisting of the redundancy (R), unique contributions from gene Y (Y) and gene Z (Z), and the synergistic contribution (S). The vertical dashed gray lines in the one-edge plot indicate the unique PID values that are used as the basis for our inference algorithm. Here, all regulatory interactions are assumed to be activating, the additional stimulating ligand targeted gene X, and the values indicated are the mean PID values calculated from five sets of simulations (with different randomly sampled initial conditions); results obtained with models that include both activating and inhibitory regulation are shown in Figure S1. (B) Visualizations of the directed 50-node networks, produced by GeneNetWeaver (Schaffter et al., 2011); node degree distributions for the 50- and 100-node GeneNetWeaver networks used in this study are shown in Figure S2. (C) Mean PID profiles for gene triplets in the 50-gene S. cerevisiae 1 network. Every triplet of nodes (genes) in the network was assigned to one of six possible classes (based on the known connectivity of genes, as indicated in network diagrams above each plot). Each line graph shows the mean PID values calculated across triplets with the same topology, with the horizontal axis labels indicating the PID contribution.

Figure 3

Demonstration of Unique Information in In Silico Networks (A) PID values are estimated from data simulated from a 10-gene in silico network (top) using GeneNetWeaver (Schaffter et al., 2011). Each line graph shows PID values estimated using genes 1 and 7 as the sources, and each of the remaining genes in turn as the target (graph titles, GX, indicate the target gene). Four PID values are given in each graph — the redundancy (R), the unique information between gene 1 and the target (G1), the unique information between gene 7 and the target (G7), and the synergy (S). The mutual information between two genes is the sum of their unique information and the redundancy (Equation 8). The ratio of the unique information to the mutual information tends to be higher between pairs of connected genes (dashed vertical lines indicate the unique contributions for connected genes). (B) Ratio of unique information to mutual information in triples with the two most common topologies, within the 50-gene network, S. cerevisiae 1. For each gene triple of the unconnected and one-edge topologies (see Table S1 for topology frequencies), we calculate the unique information values between each pair of genes and their mutual information. The ratios of unique information to mutual information are higher in general for the connected pair; the same pattern was observed in all networks. The boxes show the quartiles (the median is indicated in red) and the whiskers extend to (Q1 – 1.5 × IQR) and (Q3 + 1.5 × IQR) or to the minimum or maximum data point if these are closer.

Figure 4

The PIDC Inference Algorithm (A) Schema of the PIDC inference algorithm. PID values are estimated for every gene triplet (with each gene treated as the target gene in turn), and from these the PUC, u_X,Y, is estimated for every pair of genes. For each gene, X, an empirical distribution, f_X(u), is estimated from its PUC scores with all other genes. The confidence of an edge between a pair of genes depends on the corresponding cumulative distribution functions, F_X(u), for each gene within the pair (i.e., the blue shaded areas); these confidence scores are used to rank all possible network edges. (B) Example empirical distributions of PUC scores by gene. Gamma distributions were fitted to the PUC scores (Equation 10) for each gene in a 50-node in silico network (for each gene X, a PUC score, u_X,Y, is obtained for that gene paired with each other gene Y in the network). Due to the variability of these distributions, using a universal threshold for inferring edges is problematic, thus we use the cumulative probability distributions for each gene to obtain a final confidence score for network edges. (Colors are to aid distinguishing the distributions.)

Figure 5

Performance Comparison of Information Theory Network Inference Algorithms (A) AUPR is calculated for several algorithms applied to ten in silico datasets, generated from five 50-gene and five 100-gene networks, using Bayesian blocks discretization and the maximum likelihood estimator. The PID-based algorithm PIDC and the raw PUC score perform well in general, as does CLR. All algorithms perform better with larger datasets, but for the larger networks, this improvement is more marked in the algorithms that consider higher-order information or network context, suggesting that these are important principles for inferring networks from single-cell data. (B) Dropouts are simulated from the medium-sized 100-node datasets: the lowest 20% (low rate) or 50% (high rate) of expression values for each gene each have a 50% probability of being set to 0. Relative performance of the algorithms is the same in the presence of dropouts, although the performance of all algorithms deteriorates with a higher proportion of dropouts. (C) Violin plots of AUPR scores for all algorithms from all datasets demonstrate the value of higher-order information (PUC improves on MI; PIDC improves on CLR) and of network context (CLR improves on MI; PIDC improves on PUC); x indicates 50-gene network; ^∗ indicates 100-gene network; + indicates 100-gene network with dropouts; size indicates number of cells in the dataset. All algorithms used are described in Box 2 and Results, with MI indicating the use of mutual information scores alone to rank edges (i.e., MI relevance network); the R package minet was used for the existing inference algorithms (Meyer et al., 2008).

Figure 6

Application of the PIDC Inference Algorithm to Experimental Datasets (A) Illustration of the relationship between the three subpopulations of MEP cells: Pre-MEP cells are enriched for erythroid/megakaryocyte progenitors but still retain some potential to differentiate into other cell types (myeloid cells); E-MEP and MK-MEP cells are strongly biased toward erythroid and megakaryocyte differentiation, respectively (for details, see Psaila et al., 2016). (B) Network inferred using the PIDC algorithm from the complete set of data from Psaila et al. (2016) (the top 2.5% of edges are shown; for clarity, only nodes connected by these edges are shown). Edge colors indicate whether these edges are also detected in networks constructed using subsets of the data (comprising data from Pre-MEP cells combined with either E-MEP or MK-MEP cells). Red edges indicate those that are present in the Pre-MEP and E-MEP network but not the Pre-MEP and MK-MEP network (i.e., erythroid specific), while blue edges indicate the reverse scenario (i.e., megakaryocytic specific); edges present in both networks, or only in the network constructed using all the data, are shown in gray. (C) Comparison of the networks inferred from the data from Psaila et al. (2016) using our PIDC algorithm or MI relevance networks (in both cases, the top 2.5% of edges are shown). (D) Percentage of nodes that are connected in networks inferred from the Psaila et al. (2016) data as the threshold for edge inclusion is varied (from 0% to 100% of possible edges, according to their rank). These results show that networks inferred using PIDC (black) tend to be better connected than networks inferred using MI (red); i.e., MI networks show more skewed degree distributions. (E) PIDC interaction network inferred using early embryonic development data (oocyte to E4.25 blastocyst stages); see Guo et al. (2010). Graph edges indicate the top 5% of putative interactions detected using the PIDC algorithm on the complete dataset. Networks are also inferred using two overlapping subsets of the data: an early subset that includes all cells collected from oocyte up to 32-cell E3.5 blastocyst stages, and a late subset including cells collected from 16-cell morula to 64-cell E4.25 blastocyst stages. Edge colors indicate temporal dependencies of the identified relationships; red indicates an edge ranks in the top 5% of edges in the early network but not the late network, blue indicates the converse (late but not early), while gray indicates relationships without specific temporal dependencies (i.e., only present in the network constructed using the complete dataset, or present in both the early/late networks). (F) Comparison of networks inferred using hematopoietic development data in Moignard et al. (2015). The authors used single-cell expression data for 20 transcription factors to infer a Boolean network model of blood development; we show a simplified representation of their model, where nodes (genes) are linked by an edge if those genes are either directly linked, or linked via one Boolean operation or set of update rules (genes linked in their model via a chain of multiple sets of update rules are not connected here). We used our PIDC algorithm or MI alone to infer networks of putative interactions between genes using these same data and compare the edges identified in each of these three networks (numbers of shared edges are indicated by the Venn diagram).

Figure 7

Speed of MI and PIDC Calculations Using the Julia Programming Language (A) Comparison of the times taken to calculate a matrix of pairwise MI values using the R package minet (Meyer et al., 2008) and our Julia package InformationMeasures.jl. Input data were simulated expression values for up to 1,000 genes, with 700 values per gene (equivalent to 700 cells, the same as our medium-sized in silico dataset). Data were discretized using the uniform width algorithm, because Bayesian blocks was not available via minet, and times were measured using inbuilt functions in R and Julia. (B) Times taken to calculate the MI matrix in the Julia programming language with larger numbers of cells for networks of different sizes. Data were discretized using the recommended Bayesian blocks algorithm, which has a much greater complexity than the uniform width algorithm, but produces better estimates. (C) Time taken to infer networks of varying sizes using the PIDC algorithm implemented in Julia. Networks were inferred for simulated datasets of up to 1,000 genes, with 700 expression values per gene.