Parameter estimation in FACS-seq enables high-throughput characterization of phenotypic heterogeneity

Phenotypic heterogeneity is a most fascinating property of a population of cells, which shows the differences among individuals even with the same genetic background and extracellular environmental conditions. However, the lack of high-throughput analysis of phenotypic diversity has limited our research progress. To deal with it, we constructed a novel parameter estimation method in FACS-seq, a commonly used experimental framework, to achieve simultaneous characterization of thousands of variants in a library. We further demonstrated the model's ability in estimating the expression properties of each variant, which we believe can help to decipher the mechanisms of phenotypic heterogeneity.


Introduction
Cells are sophisticated instruments driven by the central dogma that build varieties of lives. Although traits of a cell seem simply determined by genes and environment, phenotypic diversity was found among genetically identical cells even under the same environmental condition [1]. This phenomenon, referred to as phenotypic heterogeneity, has been a major research hotspot in quantitative biology [2] that reveals non-deterministic [1] and non-linear [2] [3] properties of intracellular processes.
Hitherto, there are mainly three methods to quantitatively characterize phenotypic heterogeneity: time-lapse microscopy [4], flow cytometry [5] [6] and single-cell omics [7][8] [9]. However, each of them can only assay one type of cell at a time, which is time-consuming and labor-intensive when testing different variants. On the other hand, since the phenotypic heterogeneity itself served as a heritable trait [10], obtaining the comprehensive landscape between genetic variation and corresponding expression noise can help to deepen our understanding of the causal relationship between them. Thus, there is a critical need to design a novel method that enables high-throughput characterization of phenotypic diversity.
To address this issue, we focused on a recently developed powerful experimental framework, fluorescence-activated cell sorting and high-throughput DNA sequencing [11][12] (FACS-seq, Figure 1A). This approach enables us to associate genetic variation with the expression of the fluorescent protein, thereby allowing quantification of the expression levels of thousands of variants in parallel. However, it is generally hard to figure out the phenotypic heterogeneity of each variant. We thought this issue may root in the underutilization of data, as one conventionally less explored phenomenon is that the phenotypic heterogeneity can affect the fraction of cells across bins when sorted based on their fluorescence intensity through FACS [13] ( Figure 2B). Besides, sequencing data ( Figure 1B) and the fluorescence intensity distribution ( Figure 1C) of the whole library were neglected in previous works, which may lead to the loss of important information. Therefore, by taking these factors into consideration, we may have chances to derive the detailed properties of gene expression.
In this paper, based on the empirical discovery that gene expression follows a log-normal distribution, we present a novel parameter estimation method that enables quantification of phenotypic heterogeneity in a massively parallel manner using FACS-seq data. Our model contains two parts, in the first part we apply a maximum likelihood estimation (MLE) to get the parameter estimates for frequency distribution; while for the second part, a parametric generative adversarial network [14] (GAN) is ultilized to fit the overall fluorescence distribution. These two parts are assembled into an artificial neural network where parameters can be updated through back-propagation. As a result, this model achieves remarkable performances in both simulation and actual FACS-seq data, which paves the way for the study on sequence-heterogeneity association.

General framework of FACS-seq
FACS-seq is a method that combines two powerful high-throughput experimental frameworks together ( Figure 1A), where the FACS part can quickly examine the fluorescence intensity (x, log-scaled) of millions of cells as well as separate them based on their intensity value. Hence, an overall fluorescence distribution (p(x), Figure 1C) can be measured through it. Besides, based on the customized boundaries, cells can be sorted into several bins. Suppose in an experiment, a library with n mutants are sorted into K bins with boundaries where a cell with fluorescence intensity of x will be screened into bin k if b k−1 < x ≤ b k . The NGS part enables quantification of the abundance of each variant in a population, thus the proportion of arbitary variant i (π i , Figure 1B) in the library can be obtained. Moreover, based on the overall distribution and the sequencing data of each bin, we can further derive the frequency distribution of arbitary variant i in bin k (denoted by f ik , see Results), which satisfies (1)

Gene expression follows a log-normal distribution
The key point of describing the phenotypic heterogeneity is to find an effective statistical measurement to characterize the fluctuation in gene expression, where the coefficient of variance is commonly used in empirical studies. To derive it, the knowledge of gene expression distribution is helpful that enables precise estimation of corresponding parameters. Fortunately, plenty of works have shown that the gene expression follows a log-normal distribution [15] (Figure 2A), which can be represented using two parameters, mean and standard deviation. Suppose a library contains n variants, the log-scaled expression level of variant i (denoted by x i ) would follow a normal distribution with parameter set (µ i , σ 2 i ):

Model architecture
According to the above descriptions, the probability of sorting variant i into bin k is While the observed fraction of variant i in bin k is Therefore, a straightforward idea is to apply MLE to get corresponding parameter estimates, which is equivalent to minimize the following cross entropy loss function, whose derivatives with respect to parameters can be calculated using error function: A model structure can be simply derived based on the above loss function ( Figure 3A, red rectangle). However, the vanishing gradient problem was found when applying this model to estimate parameters in computer simulation, resulting in poor fitting for several variants. To deal with it, more detailed information should be taken into consideration to derive more accurate estimates. For this purpose, we focused on the overall fluorescence intensity distribution, which should be a mixture of Gaussians: Hence, the data generating process can be simulated using a reparameterized [16] generator G ( Figure 3A, blue rectangle; Figure 3B), which specifically includes three steps: (1) sample an arbitaty variant i from the categorical distribution P (i|π) with probability π i , which can be represented by P (z = i|π) = π i ; (2) sample a random variable u from a standard normal distribution, u ∼ N (u|0, 1); (3) obtain the sample data of log-scaled fluorescence intensity by µ i + σ i u. Moreover, to make the generator approximate to the true data distribution, a neural network-based discriminator D was applied to determine whether the data is real or fake. In other words, a two-players game is played with value function V (G, D): Combine the above two parts togather, we can conclude the whole algorithm as below 2.3.
We then fed the above information exclude µ and σ into our model. As a result, the agreement between the estimates and real parameters is very good for both mean ( Figure Figure 4A) and standard deviation ( Figure Figure 4B).

Case study: FACS-seq-based profiling of tnaC-mediated tryptophan-dependent gene expression systems
As a molecular sensor, tnaC encodes a 24-residue leader peptide that responds to intracellular tryptophan as well as regulates the biosynthesis of indole. Wang et al. [13] constructed a comprehensive codon-level mutagenesis library of tnaC based on a two-color reporter system in E. coli, where the egfp is under the control of tnaC as the sensor true } from real fluorescence intensity distribution. • Update the discriminator by ascending its stochastic gradient: end for • Sample minibatch of m samples {x (1) f ake , x (2) f ake , ..., x (m) f ake } from G(x|π, µ, σ T σ). • Update the generator by ascending its stochastic gradient: end for The agreement between estimates and real parameters is very good for both (A) mean (r 2 = 0.9986, n = 1, 500) and (B) standard deviation (r 2 = 0.9862, n = 1, 500). (C) The generated distribution is approx to the real data distribution.
response reporter, and mcherry is constitutively expressed to normalize cell-to-cell variation. Cells treated with a particular concentration of ligand were sorted into six bins according to their responses, the read number of each bin was quantified through NGS. As a proof of the concept, we tested our model on their dataset.
NGS clean data were obtained from Bioproject PRJNA503322 of NCBI Sequence Read Archive. Pairs of paired-end data were merged by FLASH script and those reads without detected pairs were removed. The regular expression pattern of 'TAAGTGCATT[ATCG]{70,80}TTTGCCCTTC' in the sequencing reads (and the reverse complementary sequence) was searched, and those carrying mutations within the upstream (TAAGTGCATT) or downstream (TTTGCCCTTC) flanking regions (10 bp each) were removed, the read number of each distinct sequence was counted.
For the two-color report system, cells were sorted based on the ratio between the log-scaled fluorescence intensity of eGFP and mCherry, Suppose a set of boundaries {eGF P = a · mCherry + b j |j = 1, ..., K − 1} were used to sort the cells into K bins, we can define x = eGF P/mCherry a to represent the expression of each cell, thus the boundaries corresponding to x should be b = {b 0 = −∞, b 1 , ..., b K−1 , b K = ∞}. As for the NGS data of each bin, the read counts were first normalized by a factor C k /R k , where C k is the total number of cells sorted in bin k and R k is the total read counts in bin k , these values were then normalized across all bins to derive the frequency distribution of each variant.
After the above pretreatment, we got three datasets corresponding to three different concentrations of ligand supplementation (0, 100, 500 µM ). We then performed our method on these datasets, as a result, the correlations between replicates proved the ability of our model that can capture the properties of gene expression of each individual in a pooled library ( Figure 5 A-F). However, the model performances in predicting the standard deviation were not so good compared to simulation results, which we thought were due to the observational errors during experiments. These errors may include polymerase preference, the random fluctuation in testing fluorescence intensity, the misclassification when sorting cells, etc. We further discovered that those errors are more likely to influence the variants with small mixing coeffecient (5 D-F), this phenomenon is obvious since few measurements of them would exacerbate the statistical uncertainty.

Conclusion
In this paper, we proposed a novel method that can precisely estimate the parameters of gene expression distribution based on the data derived from FACS-seq experiments. This model makes full use of all observation data as well as the assumption that gene expression follows a log-normal distribution, which to the author's best knowledge, is the first effective model that enables high-throughput quantification of phenotypic heterogeneity. We believe this model can help to facilitate the research progress of corresponding fields.
In this model, we only focus on the ideal situation because we don't have any prior knowledge of these observational errors. A more robust model can be constructed if we can measure the distribution of these errors, or we can reduce them in our experiment to improve data quality to achieve better performances. To derive more reliable parameter estimates, we recommend ensuring the uniformity of the library as those results associated with small proportions are more likely to be affected by the experimental noise. Besides, more bins are recommended when sorting cells as they can produce more information.

Data and code availability
The code and data related to this work can be accessed via Github (https://github.com/fenghuibao/Deep_FS.git).