Mathematical modelling of transcriptional heterogeneity identifies novel markers and subpopulations in complex tissues — TEI-tagged view

Niya WangpersName, Eric P. HoffmanpersName, Lulu ChenpersName, Li ChenpersName, Zhen ZhangpersName, Chunyu LiupersName, Guoqiang YupersName, David M. HerringtonpersName, Robert ClarkepersName, Yue WangpersName (2016date). “Mathematical modelling of transcriptional heterogeneity identifies novel markers and subpopulations in complex tissues.” Scientific Reportstitle·j.

Object of study — tissue heterogeneityterm·dom in mixed gene expression data

Research question — how to identify subpopulationterm·key marker genesterm·dom without a priori information

Contribution — CAMterm·mth, a fully unsupervisedterm·key method with proven theorems for finding marker genesterm·dom and subpopulationsterm·key

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Mathematical modelling of transcriptional heterogeneity identifies novel markers and subpopulationsterm·key in complex tissues Abstract Tissue heterogeneityterm·dom is both a major confounding factor and an underexploited information source. While a handful of reports have demonstrated the potential of supervised computational methods to deconvolute tissue heterogeneityterm·dom, these approaches require a priori information on the marker genesterm·dom or composition of known subpopulationsterm·key. To address the critical problem of the absence of validated marker genesterm·dom for many (including novel) subpopulationsterm·key, we describe convex analysis of mixturesterm·mth (CAMterm·mth), a fully unsupervisedterm·key in silico method, for identifying subpopulationterm·key marker genesterm·dom directly from the original mixed gene expressions in scatter space that can improve molecular analyses in many biological contexts. Validated with predesigned mixtures, CAMterm·mth on the gene expression data from peripheral leukocytes, brain tissue, and yeast cell cycle, revealed novel marker genesterm·dom that were otherwise undetectable using existing methods. Importantly, CAMterm·mth requires no a priori information on the number, identity, or composition of the subpopulationsterm·key present in mixed samples, and does not require the presence of pure subpopulationsterm·key in sample space. This advantage is significant in that CAMterm·mth can achieve all of its goals using only a small number of heterogeneous samples, and is more powerful to distinguish between phenotypically similar subpopulationsterm·key.

Tissue heterogeneityterm·dom, arising from multiple subpopulationsterm·key within a sample, is both a major confounding factor in studying individual subpopulationsterm·key and an underexploited information source for characterizing complex tissues 1 2 . Since the interactions among subpopulationsterm·key are fundamental to both normal development and disease progression, molecular analysis of subpopulationsterm·key in their native microenvironment provides the most biologically relevant picture of the in vivo state 3 4 . Complex tissues can be characterized by the identity, composition, and expression profile of possibly unknown subpopulationsterm·key 5 , where subpopulationsterm·key are often defined by marker genesterm·dom (genes whose expressions are exclusively enriched in a particular subpopulationterm·key 6 7 , Fig. 1a ). Current global profiling methods can neither identify differentially expressed genes among different subpopulationsterm·key, nor distinguish among the contributions of different subpopulationsterm·key to a globally measured gene expression profile 1 5 . Thus, it is generally impossible to tell whether expression change reflects a change in subpopulationterm·key composition, a change in subpopulation-specific expression, or both.

An experimental solution to mitigate tissue heterogeneityterm·dom is to isolate subpopulationsterm·key before molecular profiling by supervised cell sorting or tissue microdissection 1 8 . However, these methods are biased, costly, inapplicable to previously-assayed samples, and may alter the expression values 5 6 , restricted to cell-types that can be separated from the others. While some reports have demonstrated the potential of computational methods to resolve tissue heterogeneityterm·dom, a priori information on the composition 2 5 9 or signatures 6 10 11 12 of the subpopulationsterm·key believed to be present is almost exclusively required. Acquiring these prior information relies on experimental solutions and key limitations remain unsolved. Such supervised methods consequently have difficulty detecting subpopulationsterm·key that are subtle, condition-specific (molecular signatures and cell function are changed but not cell appearance), or previously unknown 3 13 .

To address the critical problem of the absence of validated marker genesterm·dom for many (including novel) subpopulationsterm·key, we developed a fully unsupervisedterm·key computational method (convex analysis of mixturesterm·mth – CAMterm·mth) that can identify these marker genesterm·dom directly from the original mixed expressions in scatter space-a nontrivial task. CAMterm·mth requires no prior information on the number, identity, or composition of the subpopulationsterm·key present in mixed samples 12 , and does not require the presence of pure subpopulationsterm·key in sample space 14 15 . Fundamental to the success of our approach is the newly-proven mathematical theorems, showing that the scatter simplexterm·dom of mixed expressions is a rotated and compressed version of the scatter simplexterm·dom of pure expressions, where the marker genesterm·dom are located at each vertex ( Fig. 1b ). CAMterm·mth works by geometrically identifying the vertices (and their resident genes) of the scatter simplexterm·dom of globally measured expressions (Methods).

Tissue samples to be analyzed by CAMterm·mth contain unknown numbers and varying proportions of distinct subpopulationsterm·key. Expression of a given gene in a specific subpopulationterm·key is modeled as being linearly proportional to the abundance of that subpopulationterm·key 5 6 (without log transformation 16 , Fig. 1c ). Because many genes can be co-expressed across different subpopulationsterm·key, CAMterm·mth instead identifies the marker genesterm·dom by detecting the simplex vertices of mixed expression data. The minimum description lengthterm·mth (MDLterm·mth) criterion determines the number of subpopulationsterm·key present 17 (Methods).

Results Validation of CAMterm·mth on real benchmark dataset ( GSE19830 ) Unsupervisedterm·key identification of marker genesterm·dom from mixed expression data allows us to acquire the relative expression levels of those genes ( Fig. 1a ). The average of sum-normalized marker geneterm·dom expressions produces subpopulationterm·key proportions (Methods and Fig. 1b ). Using predesigned RNA mixing experimental data acquired from biological mixtures of pure gene expressions (brain, liver, lung) 5 6 , we showed that CAMterm·mth identified the marker genesterm·dom that define each of the multiple subpopulationsterm·key ( Fig. 2a,b ) and estimated the proportions of these subpopulationsterm·key in the mixed samples ( Fig. 2c ) and their respective expression profiles ( Fig. 2d ). Since the presence of marker genesterm·dom is both a sufficient and necessary condition for deconvolutionterm·mth (Methods), these results (validated by the ground truth) confirm the existence of marker genesterm·dom and CAMterm·mth’s ability to detect these genes blindly and correctly ( Fig.2b ). Moreover, CAMterm·mth enabled detection of condition-specific marker genesterm·dom across sample groups (for example, disease versus control). Thus, novel marker genesterm·dom for a subpopulationterm·key in a given context can be determined, despite an expected change in that subpopulationterm·key’s relative abundance and/or state.

We then validated and compared the marker genesterm·dom identified by CAMterm·mth with the results of a standard differential analysis using mixed expressions, in a situation in which all ground truths are known. We reconstituted mixed expressions by multiplying the measured pure expressions (brain, liver, lung) by the proportions of subpopulationsterm·key in a given heterogeneous sample. Such mixtures mimic biological samples with varying relative abundances of different subpopulationsterm·key. Experimental results on three different mixture settings ( Fig. 3a–c corresponding to 30 ° , 45 ° and 60 ° simplex rotations) show that CAMterm·mth identified a significant majority (>90%) of ground truth marker genesterm·dom that were largely missed by standard differential analysis ( Fig. 3a–c ). Moreover, the overall constituent proportions estimated by CAMterm·mth correlate well with the ground truth ( r > 0.95; Supplementary Table S1 ). We also analyzed the same datasets with supervised methods and showed their instability in the presence of uncertainly measured cell proportions or signatures 5 6 , referring to the significantly degraded performance when the required supervising information is inaccurate or unreliable ( Supplementary Figs. S1 and S2 , Supplementary Discussion ).

Validation of CAMterm·mth on real benchmark dataset ( GSE11058 ) In most physiologically relevant settings, few universal marker genesterm·dom are available for all subpopulationsterm·key in a tissue sample, particularly under different conditions. Furthermore, no ‘gold-standard’ marker geneterm·dom list exists to accurately define many of the known subpopulationsterm·key that may be present in a sample. To test CAMterm·mth’s ability to detect novel marker genesterm·dom in complex tissues, we applied CAMterm·mth to the benchmark human blood gene expression data 11 . In this dataset, the constituent subpopulationsterm·key are four immune cells in blood that are phenotypically very similar to each other: Raji (B-cell), IM-9 (B-cell), Jurkat (T-cell), and THP-1 (monocyte). The expression profiles of the purified subpopulationsterm·key were not used in any step of CAMterm·mth but served as the ground truth for performance assessment.

Without using any a priori information, the minimum value of the MDLterm·mth curve correctly indicated the presence of four subpopulationsterm·key in the mixtures ( Fig. 4a ); CAMterm·mth blindly identified 301 marker genesterm·dom ( Fig. 4b ) and accurately estimated both the constituent proportions ( r = 0.96; Fig. 5 , Supplementary Table S2 ) and the subpopulation-specific expression profiles ( r = 0.95~0.98; Supplementary Fig. S3 ). CAMterm·mth even slightly outperformed the supervised method that used the known reference marker geneterm·dom signatures ( Supplementary Table S3 ). These results further support the existence of novel marker genesterm·dom and CAMterm·mth’s ability to detect them both blindly and effectively. Note that the correlation between the estimated and measured pure expression profiles is assessed over marker genesterm·dom rather than all genes since the high correlation obtained using all genes would be misleading due to the presence of many co-expressed genesterm·dom ( Supplementary Table S4 ). We next assessed the biological plausibility of the 301 blindly detected marker genesterm·dom. Analysis with Ingenuity Pathway Analysis (IPA) shows that these novel marker genesterm·dom are enriched in functions associated with each of the four subpopulationsterm·key ( Supplementary Fig. S4 , Supplementary Discussion ).

Application of CAMterm·mth on real benchmark dataset (yeast cell cycle) Interactions among cells are fundamental to many biological processes including repopulation dynamics 4 . As a more challenging case involving only a single cell type, we studied yeast grown under varying conditions to reinterpret cell cycle dynamics on the basis of their gene expression changes over time 10 18 . The mixed expressions were acquired from synchronously growing yeast cells in different phases of the cell cycle, where cells exist at different phases even within a given cell type, presenting an additional layer of complexity to identify differentiation marker genesterm·dom. In marked contrast to supervised methods that require a priori information, CAMterm·mth automatically detected phase-specific marker genesterm·dom of yeast cells and uncovered the subpopulationterm·key dynamics (cell cycle phase distributions) ( Fig. 6 ). These results correlate well with the patterns reported by a supervised method in the same dataset 10 . We also assessed these results by comparing with the phase-specific patterns generated by the marker genesterm·dom of Spellman et al. 18 and obtained high accordance ( Supplementary Fig. S5 , Supplementary Discussion ). Gene enrichment analysis with DAVID confirmed the biological plausibility of the 187 phase-specific marker genesterm·dom blindly identified by CAMterm·mth 19 ( Supplementary Table S5 , Supplementary Discussion ).

Discussion Most of the large datasets available for study in many diseases are based on the use of complex tissue samples with unknown subpopulationsterm·key and where no supervising information is available 12 20 . Currently, CAMterm·mth provides the only validated tool to identify marker genesterm·dom and to resolve tissue heterogeneityterm·dom in these public datasets ( Fig. 7 ). The significant advantage of our strategy is its ability to detect novel marker genesterm·dom in an unsupervisedterm·key way using the data obtained from complex tissue samples. Importantly, CAMterm·mth operates in the scatter space of the mixtures, and hence does not require the presence of pure-subpopulation samples (phenotypic archetypes) 14 15 . This advantage is significant in that CAMterm·mth can achieve all of its goals using only a small number of heterogeneous samples (when the number of heterogeneous samples equals or is greater than the expected number of distinctive subpopulationsterm·key), and is more powerful to distinguish between phenotypically similar subpopulationsterm·key (leveraging the mixing diversity across heterogeneous samples). Supported by the MDLterm·mth model selection, the benefits of CAMterm·mth method include its wide and unbiased applicability, and exciting potential for characterizing subpopulationsterm·key in tissue remodeling ( Supplementary Discussion ). While the principal application here involves gene expression data, our methodology can be readily applied to other types of quantitative measurements 21 . We have applied CAMterm·mth to unmixing other types of data such as medical imaging and remote sensing 22 23 24 25 26 and obtained very promising results.

The CAMterm·mth method is principally applicable to the situation when the number of available samples is less than the number of subpopulationsterm·key, i.e. , an underdetermined problem, albeit with some flexibility on model assumptions 27 . However, in such situation, while the CAMterm·mth solution will still be able to identify the mixing matrix (including the most appropriate number of sources), the subpopulation-specific gene expression profiles cannot be accurately estimated 27 .

Supervised methods are poorly suited to finding subpopulationsterm·key that are subtle, condition-specific (their distinctive signatures are changed when the cells are present in different microenvironments or tissues), or previously unknown. Unsupervisedterm·key approaches have not previously been reported, and yet these are necessary and urgently needed for appropriate data analysis. To provide some perspective, adequate information needed to guide supervised analyses is lacking in most of the datasets within large international studies and databases such as TCGA, METABRIC, GEO, ONCOMINE and many others.

Alternative approach by Schwartz et al. 14 and by Hart et al. 15 , which may be the closest to our approach, was not designed to address the same problems as ours and it is not an unsupervisedterm·key approach. This method operates in sample space, requires the presence of pure subpopulationterm·key samples (phenotypic archetypes), and uses a supervised/subjective scheme to determine the number of constituent subpopulationsterm·key 15 . The differences (and consequent limitations) are highly significant. Regarding the key difference between the CAMterm·mth here and alternative methods, we should emphasize that our CAMterm·mth method operates in scatter space and requires only a limited number of subpopulation-specific markers 23 24 . In scatter-space the power of detecting simplex vertices solely depends on the non-singularity (diversity) of the mixing proportions yet is independent of phenotypic similarity among subpopulationsterm·key. In contrast, in the sample-space approaches, pure subpopulationterm·key samples are required for each of the subpopulationsterm·key. While as assumed by Schwartz et al. that pure subpopulationterm·key samples may be available in the limit of uniform, dense sampling 14 , this assumption is often invalid when sample size is relatively small, and often close phenotypic similarity between subpopulationsterm·key will be a problem in identifying the vertices in sample space. These issues could greatly limit the applicability of the sample-space based unsupervisedterm·key methods.

In relation to previous work, while this is the first report on CAMterm·mth methodology with specific designs and applications to identifying novel markers and subpopulationsterm·key in complex tissues using gene expression data, the concept of exploiting the convex analysis for general unsupervisedterm·key deconvolutionterm·mth is shared by prior publications 25 . For example, sample space or compartment model based method has been applied to analyzing remote sensing or medical imaging data 22 24 26 28 , facilitated by the open-source software tool 23 29 . In contrast, the more complete mathematical derivations and CAMterm·mth algorithm, e.g. , the MDLterm·mth formulation, marker geneterm·dom concept, and theorem proofs, are specifically derived for analyzing gene expression data in biological tissue samples.

Methods Latent variable model on mixed gene expressions in heterogeneous samples Consider gene expression measured from a sample composed of K subpopulationsterm·key. We assume that the measured expression level x is the weighted sum of each subpopulationterm·key’s expression, where the contribution from a single subpopulationterm·key is proportional to the abundance and specific expression of that subpopulationterm·key 2 5 6 10 16 . The measured expression level thus is ( Fig. 1c ) where s k ( i ) is the expression level of gene i in subpopulationterm·key k , x j ( i ) is the expression level of gene i in heterogeneous sample j , and a jk is the proportion of subpopulationterm·key k in heterogeneous sample j . We further assume that gene expression values are non-negative (before log-transformation 13 16 ) and adopt the definition of subpopulation-specific marker genesterm·dom as those genes whose expression values are exclusively enriched in a particular subpopulationterm·key 6 10 11 ( Fig. 1a ). Thus, the specific expression of a marker geneterm·dom (MG) in subpopulationterm·key k * is When marker genesterm·dom are known for each subpopulationterm·key, we can use the expression values of marker genesterm·dom to deconvolute mixed expression profiles 6 10 11 . When no such prior knowledge is available ( i.e ., none of K , a jk and s k ( i MG ) is known a priori ), solving latent variable model (eq. (1) ) is essentially a blind source separation problem 22 25 , where accurate identification of subpopulation-specific marker genesterm·dom is a critical but nontrivial task 6 10 20 .

Our formulation dissects complex transcriptional heterogeneity into combinations of distinct subpopulationsterm·key, leveraging the advantages of both tissue-wide and single-cell approaches 14 30 . Specifically, discerning differences among single cells can gain valuable information about inter-cellular heterogeneity but allow only a few markers per cell and is prone to cell-cycle confounders; while tissue-wide measures provide a detailed picture of averaged population state but at the cost of losing information about inter-subpopulation heterogeneity.

Parallelism between latent variable model and the theory of convex sets Consider a set of J (≥ K ) heterogeneous samples of varying composition of unknown subpopulationsterm·key. Applying a sum-based standardization to gene expression values x j ( i ) across samples and using vector-matrix notation, we can re-express equation (1) as where x ( i ) and a k are the vector notations (over samples) of mixed expression values and subpopulationterm·key proportions, respectively. Since s k ( i ) is non-negative and standardized, as a non-negative linear combination of , the set of gene expression vectors x ( i ) forms a subset of the convex set uniquely defined by the set of 25 26 31 ( Fig. 1b ) where N is the number of genes.

Mathematical foundation for unsupervisedterm·key identification of novel marker genesterm·dom We propose the Convex Analysis of Mixturesterm·mth (CAMterm·mth) framework to exploit the strong parallelism between a linear latent variable model (eq. (3) ) and the theory of convex sets. The novel insight is that subpopulation-specific marker genesterm·dom that define latent pure subpopulationsterm·key reside at the extremities of the scatter simplexterm·dom formed by all genes, while the interior of the simplex is occupied by co-expressed genesterm·dom (whose values are linear non-negative combinations of pure subpopulationterm·key expression values) ( Fig. 1b ). We can then identify novel marker genesterm·dom by geometrically locating the vertices of the multifaceted simplex that most tightly encloses the gene expression profiles and has the same number of subpopulationsterm·key as vertices. CAMterm·mth is supported theoretically by a well-grounded mathematical framework as summarised in the following newly proven theorems (see formal proofs in Supplementary Method ).

Lemma 1 (Scatter compression and rotation) Suppose that pure subpopulationterm·key expressions are non-negative, and x (i) = a 1 s 1 (i) + … + a k s k (i) + … + a K s K (i) where a k ’s are linearly independent and non-negative, then, the scatter simplexterm·dom of pure subpopulationterm·key expressions is compressed and rotated to form the scatter simplexterm·dom of mixed expressions whose vertices coincide with a k ’s.

The terms ‘compressed and rotated’ are used here to illustrate the geometric transformation on the simplex made by the mixing process in scatter space. Given the fact that in this case the proportions a k ’s are also non-negative (in addition to gene expressions being non-negative), as the result, Lemma 1 shows that every mixed gene expression data point x ( i ) is confined within the simplex (convex hull) defined by a k ’s. Since a k ’s are resided within the first quadrant while the simplex of pure subpopulationterm·key expressions s ( i ) originally occupy the first quadrant, we therefore use the terms ‘compressed and rotated’ to describe such transformation.

Theorem 1 (Unsupervisedterm·key identifiability) Suppose that pure subpopulationterm·key expressions are non-negative and subpopulation-specific marker genesterm·dom exist for each constituting subpopulationterm·key, and x (i) = a 1 s 1 (i) + … + a k s k (i) + … + a K s K (i) where a k ’s are linearly independent and non-negative, then, the vertices of the scatter simplexterm·dom of mixed expressions host subpopulation-specific marker genesterm·dom and coincide with a k ’s that can be readily estimated from marker geneterm·dom expression values with appropriate rescaling.

From Lemma 1 and Theorem 1, there is a feasible mathematical solution to identify subpopulation-specific marker genesterm·dom from the measured gene expression mixtures: in principle, under a noise-free scenario, we can blindly identify novel marker geneterm·dom indices by locating the vertices of the mixed expression scatter simplexterm·dom 22 23 24 . We emphasize that CAMterm·mth can distinguish between phenotypically similar subpopulationsterm·key, by working in scatter space in which the power of detecting simplex vertices depends solely on the mixture diversity (a basic requirement for any inverse problem) rather than phenotypic diversity 14 . This is significantly different from the relevant approaches (including our own other previous work) that operate in sample space and require the presence of pure subpopulationterm·key samples 14 15 25 .

Data preprocessing First, we eliminate genes whose signal intensity (vector norm) is lower than 5% (noise) or higher than 95% (outlier) of the mean value over all genes. The signals from these genes are unreliable and could have a negative impact on the subsequent analyses ( Supplementary Method ). Second, when J ≫ K , dimension reduction is performed on the raw measurements using principal component analysis, sample clustering or nonnegative matrix factorization techniques, to improve the efficiency of subsequent analyses 13 23 .

Aggregation of gene expression vectors To further reduce the impact of noise/outlier data points (gene expression vectors over samples), improve the efficiency of CAMterm·mth, and permit appropriate parameterization of the MDLterm·mth criterion to determine the number of subpopulationsterm·key, we aggregate gene vectors into representative clusters using affinity propagation clusteringterm·mth (APCterm·mth) 23 24 26 32 ( Supplementary Method ). As an initialization-free and near-global-optimum clustering method, APCterm·mth simultaneously considers all gene vectors as potential exemplars and recursively exchanges real-valued ‘messages’ between gene vectors until a high-quality set of exemplars and corresponding clusters gradually emerge. The APCterm·mth algorithm is data-driven, so the message-passing procedure may be terminated after a fixed number of iterations or after the updates stay constant for some number of iterations. In all of our experiments, we adopted a default damping factor of 0.5. The update rules are repeated iteratively and terminated when no further change occurs for about 10 iterations 23 24 . Our experience indicates that these default algorithmic parameter settings are quite suitable for obtaining good results.

Convex analysis of mixturesterm·mth (CAMterm·mth) To identify the vertices of clustered convex set X (scatter simplexterm·dom of mixed expression profiles), we performed CAMterm·mth on the obtained M cluster centers of gene vectors. We assumed K true vertices and conducted an exhaustive combinatorial search (with total combinations), based on a convex-hull-to-data fitting criterion, to identify the most probable K vertices. We used the margin-of-error to quantify the distance between and convex set X defined by , where we have if is inside X . We then selected the most probable K vertices when the corresponding sum of the margin-of-error between the convex hull and the remaining “exterior” cluster centers reaches its minimum 22 23 26 : Subsequently, we identified the indices of subpopulation-specific marker genesterm·dom based on the memberships associated with , where denote the cluster indices of the true simplex vertices, and the genes assigned to gene cluster at a vertex are declared to be marker genesterm·dom, i.e. , .

Model selection procedure One important discovery step for CAMterm·mth (as a fully unsupervisedterm·key method) is to automatically detect the number K of cell subpopulationsterm·key in the heterogeneous samples. We used MDLterm·mth, a widely-adopted and consistent information theoretic criterion 17 , to guide model selection 26 ( Supplementary Method ). We performed CAMterm·mth on several competing candidates, and selected the optimal model that assigns high probabilities to the observed data with parameters that are not too complex to encode 17 . Specifically, a model is selected with K subpopulationsterm·key by minimizing the total description code length defined by 26 where denotes the joint likelihood function of the clustered latent variable model, X M denotes the set of M gene vector cluster centers, and denotes the set of freely adjustable parameters in the clustered latent variable model 23 24 26 . Moreover, when estimating the mixing matrix ( i.e. , the column vectors a k ’s ) that is parameterized by ( K -1) J independent entries, we use some form of vector-average operation ( i.e. , ), where the scalar entry a k ( j ) is estimated from M scalar entries g m ( j ) for a given j , contributing total ( K -1) J log( M )/2 bits. Similarly, when estimating the subpopulation-specific gene expression profiles (the row vectors s k ) with total KM entries, we use a vector-average operation ( i.e. , solving linear equations), where the scalar entry s k ( m ) is estimated involving only J scalar entries g m ( j ) for a given m , contributing total KM log( J )/2 bits.

Estimation of constituent proportions and subpopulation-specific expression profiles On the basis of the expression levels of subpopulation-specific marker genesterm·dom detected by CAMterm·mth, the relative proportions of constituent subpopulationsterm·key are estimated using standardized averaging, where MG- k is the index set of marker genesterm·dom for subpopulationterm·key k ; n MG- k is the number of marker genesterm·dom for subpopulationterm·key k ; and ||.|| denotes the vector norm ( L 1 or L 2 ). The resulting are then used to deconvolute the mixed expressions into subpopulation-specific profiles by non-negative least-square regression techniques 6 10 11 23 26 .

Performance evaluation criteria We use four complementary evaluation criteria and the ground truth to assess the performance of CAMterm·mth. To assess the membership match and mismatch between the subpopulation-specific marker genesterm·dom detected from pure subpopulationterm·key expressions versus from original mixed expressions, we use Venn diagrams to visualize the outcomes and also calculate both sensitivity and specificity. To assess the accuracy of subpopulationterm·key proportion estimates , in addition to classic Pearson correlation coefficient , we also adopt the E 1 criterion given by 33 where p ij is the ij th element of the matrix . Note that E 1 is invariant to permutation or scaling, and E 1 = 0 when the estimation is perfect. Moreover, to assess the accuracy of estimated subpopulation-specific expression profiles , we calculate Pearson correlation coefficient between the estimated expression profile and ground truth over ‘subpopulation-specific marker genesterm·dom’ and ‘all genes’, respectively.

Robust detection of marker genesterm·dom on numerically mixed expressions (synthetic data) To validate CAMterm·mth’s ability to blindly detect marker genesterm·dom from mixed gene expression profiles, we numerically generated three different sets of mixtures using the profiles of pure liver, brain and lung tissues 5 , corresponding to 30 o , 45 o and 60 o simplex rotations and arbitrary compressions ( Fig. 3 ). We performed standard differential analysis (one-versus-everyone 7 ) on the gene expression profiles of the three source tissues to construct a ‘gold standard’ set of 884 marker genesterm·dom. We applied CAMterm·mth to these mixtures and compared the marker genesterm·dom blindly detected by CAMterm·mth to the gold standard set, assessed by both sensitivity and specificity. In these experiments, we set K = 3 as supported by MDLterm·mth based model selection. Using these blindly detected marker genesterm·dom, we estimated the constituent proportions in each of the mixtures and subpopulation-specific expression profiles ( Fig. 3 , Supplementary Table S2 ).

Performance comparison between CAMterm·mth and supervised methods (synthetic data) To compare the performance between unsupervisedterm·key CAMterm·mth and supervised methods, we analyzed the same benchmark datasets with supervised methods 5 6 . Since a priori information used in supervised methods is prone to various errors and is often condition-specific, we mainly assessed the reliability of these methods by numerically perturbing either the baseline constituent proportions 5 or the reference signatures 6 . We first randomly perturbed the constituent proportions with ± 0.05 ∼ 0.3 or replaced one of the mixtures with arbitrary constituent proportions, and compared the estimated subpopulation-specific profiles against ground truth using correlation coefficient over marker genesterm·dom 5 . We next randomly perturbed the reference signatures with ± 1 ∼ 50% expression variations or replaced some of the marker genesterm·dom with arbitrary ones, and compared the estimated constituent proportions against ground truth via correlation coefficient 6 . These comparison experiments were performed on a large number of replications ( Supplementary Figs. S1–S2 , Supplementary Tables S6–S7 , Supplementary Discussion ).

Accurate detection of marker genesterm·dom on biologically mixed expressions ( GSE19830 ) To further validate CAMterm·mth’s ability to blindly detect marker genesterm·dom from mixed gene expression profiles, we analyzed the benchmark real gene expression dataset GSE19830 5 . The dataset contains 3 source tissue samples and 11 biologically mixed heterogeneous samples with known mixing proportions, with each sample having 3 replicates 5 ( Supplementary Table S6 ). Using the marker genesterm·dom blindly detected by CAMterm·mth, we estimated the constituent proportions in each heterogeneous samples and subpopulation-specific expression profiles ( Fig. 2a–d , Supplementary Fig. S6 ). Besides the Venn diagrams comparing the marker genesterm·dom detected by CAMterm·mth with the gold standard ( Fig. 2b ), we also assumed the coexistence of distinct subpopulationsterm·key and subpopulation-specific marker genesterm·dom, since a distinct subpopulationterm·key would logically require a set of marker genesterm·dom that are ‘unique’ to that subpopulationterm·key. In our newly proved two theorems, the existence of subpopulation-specific marker genesterm·dom is both a sufficient and a necessary condition for a successful blind deconvolutionterm·mth. Thus, if a successful deconvolutionterm·mth of mixed expression data is achieved, the existence of subpopulation-specific marker genesterm·dom is almost certain ( Supplementary Discussion ).

Effective detection of marker genesterm·dom on biologically mixed expressions ( GSE19380 ) We then applied CAMterm·mth to the benchmark real gene expression dataset GSE19380 6 . The three rat primary subpopulationsterm·key (neuron, astrocytes, oligodendrocytes) were cultured separately, and their mRNAs were extracted and variably mixed to generate gene expression data using Affymetrix Rat 230 2.0. Without using any prior information, CAMterm·mth blindly detected 183 probe sets that are located at the scatter simplexterm·dom vertices of mixed expression profiles. These probe sets were confirmed by the classic analysis using the pure subpopulationterm·key expression profiles. We accurately estimated the subpopulationterm·key proportions with a correlation coefficient of 0.99, and the subpopulation-specific gene expression profiles with a correlation coefficient of 0.94–0.99, when compared with the ground truth. We found that while a priori marker genesterm·dom ( Supplementary Table S7 ) used by supervised method 6 were near the vertices of pure subpopulationterm·key expression scatter simplexterm·dom, some of these a priori marker genesterm·dom clearly deviated from the vertices of the mixed expression scatter simplexterm·dom.

Analyses of benchmark blood cell mixed expressions ( GSE11058 ) The dataset contains 24 samples of four subpopulationsterm·key, each with three replicated gene expression profiles of variable combinations 11 ( Supplementary Table S2 ). We averaged the three replicates on each of 12 pure and 12 mixed expression profiles, resulting in 4 pure subpopulationterm·key and 4 mixed gene expression profiles. After pre-processing and normalization by MAS 5.0, total 13,000 probe sets were retained. The 301 novel marker genesterm·dom blindly identified by CAMterm·mth from mixed expressions meet almost perfectly to the definition of marker genesterm·dom, i.e., genes whose expressions are exclusively enriched in a particular subpopulationterm·key and functionally enriched 6 7 ( Supplementary Fig. S4 ). We also compared the deconvolutionterm·mth results of CAMterm·mth and supervised methods 11 12 13 ( Supplementary Tables S3 and Table S4f ). This is particularly encouraging, since it means that novel marker genesterm·dom for phenotypically similar subpopulationsterm·key, given sufficient mixture diversity, can be detected efficiently by CAMterm·mth.

Analyses of benchmark yeast cell cycle time-course data (CDC28) Cells were collected at 17 time points taken at 10 minute intervals, covering nearly two full cell cycles of synchronized yeast cell growth under variable temperatures 10 . The gene expression profiles contain 6208 transcripts across 4 discrete phases 34 . Note that here CAMterm·mth exploited the temporal diversity to perform the unsupervisedterm·key analysis. CAMterm·mth identified total 187 phase-specific marker genesterm·dom 34 35 , among which 114 (61%) are previously reported 18 36 and 73 (39%) are considered novel, with 22 genes (21.2%) being verified cell cycle regulated genes 35 36 . Using DAVID to explore gene set functions, we found that 21 genes are enriched in G1 phase, 58 in S phase, 16 in G2/M phase, and 68 in M/G1 phase ( Supplementary Table S5 ). The repopulation dynamics curves revealed blindly by CAMterm·mth ( K = 4 suggested by MDLterm·mth, Supplementary Fig. S7 ) are highly consistent with the known cell cycle patterns 18 ( Supplementary Fig. S5 , Supplementary Discussion ). This is particularly appealing in that when a sufficient number of diversely-mixed samples is available, CAMterm·mth should uncover all subpopulationsterm·key of different cell types and identify the same cell type at different cell-cycle phases 30 .

Availability of CAMterm·mth software and supporting data A Java-R package of CAMterm·mth is available for academic and non-commercial use at: http://mloss.org/software/view/437 . In addition, public gene expression data analyzed in this paper are also available from the Gene Expression Omnibus Database under Accession Number GEO: GSE19830 , GSE19380 , GSE11058 , and CDC28.

Additional Information How to cite this article : WangpersName, N. et al. Mathematical modelling of transcriptional heterogeneity identifies novel markers and subpopulationsterm·key in complex tissues. Sci. Rep. 6 , 18909; doi: 10.1038/srep18909 (2016).

Supplementary Material