A top-down measure of gene-to-gene coordination for analyzing cell-to-cell variability

Recent technological advances, such as single-cell RNA sequencing (scRNA-seq), allow the measurement of gene expression profiles of individual cells. These expression profiles typically exhibit substantial variations even across seemingly homogeneous populations of cells. Two main different sources contribute to this measured variability: actual differences between the biological activity of the cells and technical measurement errors. Analysis of the biological variability may provide information about the underlying gene regulation of the cells, yet distinguishing it from the technical variability is a challenge. Here, we apply a recently developed computational method for measuring the global gene coordination level (GCL) to systematically study the cell-to-cell variability in numerical models of gene regulation. We simulate ‘biological variability’ by introducing heterogeneity in the underlying regulatory dynamic of different cells, while ‘technical variability’ is represented by stochastic measurement noise. We show that the GCL decreases for cohorts of cells with increased ‘biological variability’ only when it is originated from the interactions between the genes. Moreover, we find that the GCL can evaluate and compare—for cohorts with the same cell-to-cell variability—the ratio between the introduced biological and technical variability. Finally, we show that the GCL is robust against spurious correlations that originate from a small sample size or from the compositionality of the data. The presented methodology can be useful for future analysis of high-dimensional ecological and biochemical dynamics.

Biological functions within the cell are carried out by gene products, such as RNA chains and proteins 1 . The high-dimensional data of the gene expression profile is thus an elemental and useful representation of the cellular activity 2 . Until recently, only bulk RNA sequencing technologies were available to study gene expression patterns at the population level. Such measurements reflect the averaged gene expression across thousands of cells 3,4 . The recent development of single-cell RNA sequencing (scRNA-seq) technologies allow the dissection of gene expression at single-cell resolution [5][6][7][8][9] . A central observation from such single-cell measurements is that, even within populations of cells from the same tissue and of the same cell-type, the gene expression profiles substantially deviate across different individual cells.
This cell-to-cell variability may arise from two fundamentally different types of processes: (1) processes that happen while the cell is still alive and active, and (2) processes that take place while the cell content is dissected and measured. Accordingly, the cell-to-cell variability that stems from such processes can be termed 'biological' or 'technical' variability, respectively. Biological variability is the result of various sources, such as, inherent stochasticity in the biochemical process of gene expression [10][11][12] , random genetic and epigenetic mutations in different individual cells 13,14 , differences in the internal states in the cell cycle progression 15,16 , or sub-populations of cells owing to subtle environmental differences or cell-subtypes [17][18][19] . In contrast, technical variability represents statistical and measurement limitations of the single-cell sequencing procedure, e.g., sampling noise in the genetic survey and stochastic over-and under-amplification of random genes 20,21 . Naively, measuring cell-to-cell variability by evaluating the dissimilarity between the measured gene expression profiles of different cells does not discriminate between biological and technical variability. Since the real biological phenomena are manifested only in the biological variability and not in the technical noise, there is a practical need to be able to distinguish between these types of variability.
To address this challenge, several experimental and computational techniques have been developed. Some of the sources for biological variability can be controlled by dividing heterogeneous populations of cells into subpopulations which are homogeneous with respect to a particular factor [22][23][24] . For example, focusing on cells from the same cell-cycle stage or cells belonging to the same sub-type. Yet, even in such seemingly homogeneous subpopulations, biological variability still occurs from processes of stochastic nature that affect each cell differently 25 26 . A central approach considers that technical variability of individual genes follows a typical form of a variance-to-mean ratio whereas larger variability in an individual gene may indicate a biological source [27][28][29] . Other computational approaches assume that biological variability is typically characterized by interrelations between the expression of different genes, and thus, employ the measure of the correlation between pairs of genes to identify biological variability 30 . The correlation-based approaches are commonly implemented in a 'bottom-up' fashion, e.g., by calculating co-expression matrices 31 . However, without a priori knowledge of the intricate map of gene-to-gene regulatory interactions, it is extremely difficult to accurately infer the interactions for a large number of genes, mainly since large calculated co-expression matrices contain a considerable amount of noise. In addition, different co-expression measures are designed to capture specific features which are not necessarily optimal for depicting all types of gene-to-gene interrelations (e.g., Pearson correlations represent only linear relationships), and are focused on pairwise interactions while an individual gene may be controlled by a combination of multiple regulators 32,33 .
Recently, a 'top-down' approach has been introduced to analyze scRNA-seq data by evaluating the global coordination level between genes (named GCL) 34 . Here, following the above-mentioned approaches that focus on the interrelations between genes to detect biological activity, we propose to use the GCL as an indication for the biological origins of the measured cell-to-cell variability. We systematically analyze synthetic data generated from mathematical models of gene regulatory dynamics, where biological variability is introduced as random variations between the generating models of different individual cells. We show that the GCL is an effective tool in the analysis of cell-to-cell variability in single-cell data. The GCL is not negligible wherever the gene regulatory models include interactions between the genes, and decreases as the amount of random variations in the dynamics increases. We also calculate the GCL for cells generated from models that involve both biological and technical variability, where the latter is introduced as independent measurement noise. We find that the GCL can distinguish between different assemblages of cells with the same measured cell-to-cell variability but a different ratio of biological-technical variability. Finally, we show that the GCL is robust against spurious correlations that originate from a small sample size or the compositionality of the data.

Methodology
Annotations used in the manuscript. In the following, we use these definitions: the gene expression of gene i in the mathematical model of gene regulatory dynamics is noted as x i . The solution of this model, which represents the steady state or the actual gene expression of gene i is noted as x * i . The measured expression, which includes both the actual expression level and measurement noise, is noted as x i .

Global coordination level (GCL).
The GCL is a "top-down" computational method to evaluate the system-wide transcriptional multivariate dependency of genes, without inferring the whole network of pairwise correlations 34 . We refer to a set of M measured cells with N genes as a matrix X N×M , where every vector column x (ν) (ν = 1 . . . M) represents the individual cell ν , and each element x (ν) i (i = 1 . . . N) represents the measured expression value of gene i in that cell. The GCL calculation is performed as follows. First, we divide the matrix X N×M into two random complementary parts A and B, where each part contains N/2 rows, i.e., each part contains N/2 gene expression values for M cells. Second, we calculate the "bias-corrected distance correlation" (bcdCorr) measure 35 on (A, B). In brief, the bcdCorr, a refined version of the "distance correlation" (dCorr) measure 36 , evaluates the level of dependence between two high-dimensional variables by testing how the distance between two samples with respect to one variable is changed compared to the distance between the same two samples with respect to other variable. Thus, the bcdCorr is a measure of the dependency level between gene-sets A and B. Finally, we repeat these two steps m times and define the GCL as the average bcdCorr, i.e., the GCL is defined as where all m divisions (A k , B k ) are independent. As the GCL typically stabilizes for m > 10 , in our analysis we choose m = 50 . For a large sample size, the empirical GCL is zero in the case of independent gene expression; whereas a significantly non-zero GCL reflects coordinated transcriptional expression, which could be interpreted as a result of underlying molecular dynamics, such as gene-to-gene regulatory interactions.

Numerical model for synthetic gene expression data.
To investigate the ability of the GCL to reveal biological and technical variability in gene expression data, we apply it to synthetic cells generated from models of gene regulatory dynamics. The model simulates 'cohorts' of gene-expression profiles with varying levels of cell-to-cell variability that originated from both differences between the actual expression profiles ('biological variability') and stochastic measurement noise ('technical variability').
We define the vector x * (ν) as the actual gene expression of an individual cell ν . The expression profile x * (ν) is modelled as the steady state of a set of coupled ordinary differential equations (ODEs), representing the gene regulatory dynamics 33,37,38 . Specifically, we use the following set of ODEs,  i,j (see Fig. 1a). In our simulation we use GRNs with self regulation and random links between the nodes, i.e., each pair of genes are connected with a constant probability in the form of an Erdős-Rényi network with an average degree equals to two. Finally, we set B = 1 , n = 1 , and the GRN weights w i,j (for existing links) are randomly selected from the uniform distribution U(0, 2).
We define the matrix X * N×M as a 'cohort' of M expression profiles x * (ν) with N genes each. For a given GRN dynamics, defined by the matrix of weights w (0) , we generate different expression profiles by performing two types of changes in the dynamics. First, a random subset of the genes in each cell are set as inoperative, i.e., their expression levels are set to zero. Specifically, in our simulations we randomly choose 5 out of 200 genes to be inoperative. Second, the GRN dynamics of each single-cell ν , w (ν) ( ν = 1 . . . M ), is generated as random variations of w (0) , where the GRNs' heterogeneity is controlled by the parameter p (see Fig. 1a). Specifically, i are generated from a normal distribution N (0, σ 2 ). To summarize, in our simulations the cell-to-cell variability of a measured cohort of cells X N×M (of a specific w 0 ) is determined by two parameters: p and σ . Thus, the 'biological variability' is generated using the parameter p which controls the heterogeneity of the underlying regulatory dynamics, and the 'technical variability' is generated using the parameter σ which controls the level of the stochastic noise.

Results
We start by demonstrating that by applying GCL on a cohort of steady-states (samples), it can capture the presence of gene-gene interactions in the underlying model. This is in contrast with the cell-to-cell variability. We compare two different models for generating cohorts of samples, X * N×M , that represent the actual gene expression of M cells with N genes, without adding measurement noise. The first model includes both self-regulation and gene-gene interaction, as detailed above in "Methodology" , while the second model has no gene-gene interactions (i.e. w (ν) i,j = 0 for any i = j in Eq. (2)). In both models, increasing the GRNs' heterogeneity level, p, increases the cell-to-cell variability (Fig. 1b,d). This is expected as the heterogeneity reduces the similarity between the equations which regulates the different cells, leading to a larger variability in the steady states. However, contrary to the variability score, the curve of the GCL score as a function of p behaves differently for these two models (Fig. 1c,e). In the first model, where the GRN dynamics contains gene-gene interactions, the GCL is significantly larger than zero for small values of p, see Fig. 1b. In addition, as the heterogeneity level increases, the gene-gene interaction are less consistent across different cells, leading to decreased GCL. In marked contrast, where the GRN dynamics does not contain gene-gene interactions, the GCL score is around zero for any value of p, see Fig. 1e. These results demonstrate that the GCL can reveal essential features of the underlying regulatory dynamics (the presence of gene-gene interactions), which are not captured by the standard measure of cell-to-cell variability.
Next, we generate and analyze data with both biological and technical variability, which are determined by the parameters p and σ , respectively. These simulations represent the measured gene expression profiles, X N×M , described above in "Methodology". We ask, given two cohorts with the same measured cell-to-cell variability, is it possible to differentiate between the one that was generated with a higher ratio of 'biological' compared to 'technical' variability? To address this question, we generated 420 cohorts, each of M = 100 cells and N = 200 genes, generated with 0.25 ≤ p ≤ 1 and 0 ≤ σ ≤ 0.5 . For each generated cohort, we calculated both the cell-tocell variability and the GCL. The heat-map in Fig. 2a shows that the cell-to-cell variability increases monotonically with each of the parameters p and σ , where different combinations of these values can lead to the same variability. Each dashed black contour line marks cohorts with equal variability, where the red line marks the cohorts with variability equals to 0.5, as a specific example. Fig. 2b shows the GCL values calculated for the same cohorts as in Fig. 2a  www.nature.com/scientificreports/ GCL calculated along the red line with respect to log(p/σ ) . Qualitatively similar behavior is also seen for cohorts where the cell-to-cell variability is equal to 0.45 and 0.55. When comparing cohorts with the same measured cell-to-cell variability, the one with a higher ratio of biological versus technical variability has a higher GCL value. These results demonstrate the inability of the cell-to-cell variability measure alone in detecting essential features of the underlying dynamics, such as distinguishing between 'technical' versus 'biological' noise. A joint analysis by both measures, i.e., cell-to-cell variability and GCL, is recommended. In addition, we compare the 'top-down' GCL and the classical 'bottom-up' co-expression matrix. The average of the gene co-expression matrix, is defined as �C� = where a matrix element C i,j represents the Spearman correlation between gene i and gene j. These two approaches were recently compared on real transcriptomic data of aging cells 34 . There was a consistent pattern of reduced GCL values in aging cells across different cell types and different organisms. In contrast, there was no clear pattern of change of the average co-expression values in old cells compared with young cells (see SI of ref. 34 ). Here, we study the effect of two typical features of real transcriptomic datasets on the ability of two approaches to reliably identify the interrelations between genes. The first feature, which is a typical scenario in currently available transcriptomic data sets, is a small sample size, i.e., the number of cells is relatively smaller compared with the number of genes and the number of possible gene-gene interactions. The second feature is the compositionality of the relative abundance of mRNAs in genomic survey data, which may lead to spurious correlations between genes [39][40][41][42] .
To study the effect of small number of cells, we first generate expression profiles X N× M with M = 150 cells and N = 200 genes, following the same procedure detailed above in "Methodology" , with p = 0.5 and σ = 0 . We then test the consistency of C and GCL values across cohorts with a decreasing number of cells. Fig. 3a shows that C becomes larger for smaller number of cells, while Fig. 3b shows that the GCL score is relatively consistent across the cohorts. Thus, we conclude that compared to the C score, the GCL is less affected by false correlations that appear due to a small number of cells. This simplified model demonstrates the disparity between the 'top-down' and the 'bottom-up' approaches. The averaged co-expression matrix, unlike the GCL, accumulates the noise from all matrix elements. This effect becomes especially pronounced in the case of a small sample size.
To study the effect of spurious correlations in compositional data, we generated cohorts of normalized 'expression profiles' with no real correlations between the genes. We generate the profiles as follows: first, we create a 'master profile' , y (0) , by generating for each 'gene' a random number from a power-law distribution with an exponent γ , i.e., p(x) ∼ x −γ . Second, we generate M = 100 profiles, where each profile i is a random number from a normal distribution with mean 1 and σ = 0.2 . Finally, we normalize each profile to 1. For each cohort of profiles, with different values of γ , we calculate both the GCL and C . Figure 3c,d show that the normalization procedure leads to spurious correlations between the genes for small values of γ , where the cells are more heterogeneous. In contrast, the GCL score is stable against this effect and correctly identifies that there is no coordination between the genes.

Discussion
To study the sources of cell-to-cell variability, we adopt the GCL method that measures the coordination between genes in a top-down approach and compared it to other classical measures. We calculate the GCL for simulated gene expression profiles, with different levels of inconsistency in the gene regulatory dynamics (as biological variability) and different strengths of measurement noise (as technical variability). We demonstrate that positive GCL values reflect the effect of interactions between the genes. We show that in the case where the variability stems only from inconsistent dynamics across the cells, the GCL decreases as the inconsistency level of gene-gene interactions increases. However, in the case of both inconsistent dynamics and measurement noise, we show that for cohorts with the same measured variability, the one with the lower GCL has lower biological variability (and higher measurement noise) compared to the other cohort.
These results may have practical applications when comparing different data-sets for studying the source of the variability. A common task in biological experiments is to compare the gene expression between two states, e.g. control vs. disease or before and after perturbation. The GCL measure allows the detection of changes in the underlying regulatory dynamics even when the mean expression values and the cell-to-cell variability do not change.
Even with the recent explosion in available transcriptomic data, we are still far away from having a sufficiently large sample size compared with the huge number of genes. Thus, in order to get a detailed description of the complex network of gene interaction in the cells, we have to be creative in our analysis. A top-down approach indeed ignores the details of the delicate interactions, but as shown here, the GCL has considerable advantages. It can help us analyze the variability source and it is robust against spurious correlations that originate from a small sample size or the compositionality of the data. We propose the GCL measure as a new tool, which we believe can provide additional insights to the classical bottom-up approach.
Practical guidelines for applying the GCL method in single-cell analysis. As discussed in this manuscript, the GCL method intends to analyze seemingly homogeneous cohorts of single-cell transcriptomes, where the measured variability stems mainly from intrinsic biological variability or technical errors. Here we suggest several pre-processing steps that are recommended to ensure that the analyzed cohorts are as homogeneous as possible, followed by a discussion of possible interpretations of the GCL outcomes.
When analyzing a set of expression profiles, the first step would be to reduce the heterogeneity by focusing on sub-populations of cells according to available metadata or specific transcriptional signatures. For example, three cell types that were isolated from a mixed population of hematopoietic stem cells using the expression of specific markers 43 , were analyzed separately by the GCL in ref. 34 . Another example for cell filtering is to reduce heterogeneity that stems from the cell cycle by selecting non-cycling cells or cells belonging to the same phase. A www.nature.com/scientificreports/ second step would be to reduce heterogeneity by performing unsupervised clustering analysis to the expression profiles and selecting only cells that belong to the same cluster or 'sub-type' . A third step is to remove outliers, i.e., cells for which their expression profile is extremely different from the average cell. A possible outlier filtering is to calculate the average profile and all distances between it and each expression profile, and then to remove cells for which the distance from the average profile is more than two standard deviations larger than the mean distance. Finally, the GCL is also susceptible to the presence of cell-pairs with very similar profiles. In real transcriptomic data, this could be due to cells that were divided just before the moment of measurement. A simple way to filter out those cell-pairs would be to calculate cell-to-cell distances between all cell pairs and remove one of these cells if the distance is exceptionally small. www.nature.com/scientificreports/ After these steps, the GCL may be interpreted as follows. When applied on a single cohort, a significant positive GCL value reflects the effect of interactions between the genes. The significance, in this case, can be evaluated by performing a Jackknife procedure on the real data and compare the results to shuffled data where the effects of interactions between genes are removed.
When two cohorts are compared, the GCL values should be investigated with regards to the measured cell-tocell variability of these cohorts. If the variability is preserved but the GCL is different, this may indicate that the underlying dynamics in the cohort with the higher GCL are more consistent across cells. For example, in ref. 34 , reduced GCL values were found to be associated with aging cells and with increased genetic mutational load, and were interpreted as random aberrations of the cellular regulatory mechanisms. However, if the variability is not the same, the GCL should be interpreted with more caution. If the lower GCL is measured in the one with the higher variability, then the GCL may provide no additional insights. This is because both increased biological variability and increased technical variability lead to lower GCL (as shown in Figs. 1 and 2 in our manuscript). But, if the GCL difference is in the opposite direction, i.e., a lower GCL is measured in the cohort with the lower variability, this may indicate a lower coordinated biological activity compared with the other cohort.

Code availability
The custom MATLAB code for computing the GCL that was used in this study is available at https:// github. com/ guy531/ gcl.