In silico microdissection of microarray data from heterogeneous cell populations

Inversion of sample heterogeneity

The first goal is to invert the mixing effect caused by sample heterogeneity. We apply the linear model developed in Materials and methods Section to the heterogeneous microarray data. The obtained results are presented below.

Because of the inherent variability of individual gene expression values, the performance of the inversion method cannot be estimated based on results for individual genes. (For illustration purposes, the results of inversion of the mixing effect for individual genes are discussed and shown later on in connection with Figures 6 and 7.) Thus, we examine the performance of our method globally, by comparing the measured and estimated expression values of all the genes simultaneously. For performance evaluation and visualization purposes, the dimensionality of the 4704-dimensional expression profiles is reduced using the standard principal component analysis (PCA). The effect of the sample heterogeneity is the same for all the genes within one array. Therefore, for each array, it is useful to combine the results over all the genes. In other words, instead of looking at individual genes, we combine the expression values of all the genes and visualize the results using the most significant principal components. For comparison purposes, we also show the samples used as a reference in the conducted microarray experiments. Since the number of measurements is far smaller than the number of genes, we use a standard approach when solving the PCA eigenvector-eigenvalue problem. Let and , i = 1, ..., K, denote the measured mixture and reference samples, respectively, and let and denote the estimated RKO and normal expression profiles. Let , where . Instead of finding the eigenvalues of the original sample covariance matrix D ^T D, we compute them for the matrix DD ^T . The eigenvalues of D ^T D and DD ^T are the same and the eigenvectors of D ^T D can be obtained from the eigenvectors of DD ^T by multiplying them by D ^T . Results of the inversion of the sample heterogeneity are shown in Figures 1 and 2. In Figure 1, all five heterogeneous samples are used to estimate the expression values of the pure colon cancer and lymphocyte samples. The two most significant PCA components of all the heterogeneous samples, reference samples, and the estimated expression profile of the pure colon cancer cells and lymphocytes are shown. Figure 1 clearly shows that the heterogeneous samples ('m1' through 'm5') are located almost on a straight line in the 2-dimensional PCA space. Furthermore, the line on which the heterogeneous samples are lying is parallel to the first principal component, suggesting that the most significant variation in the data is due to the linear mixing effect. The estimated expression profile of the pure colon cancer cells and lymphocytes are close to samples number #1 and #5, respectively, indicating that the inversion of the mixing phenomenon produces reasonable results.

The results are more easily appreciated when only the most significant PCA component is shown. As discussed above, the variation in the most significant PCA component is due to the mixing effect. The results in Figure 2 (a) are as in Figure 1, but now shown in 1-dimension in order to facilitate the interpretation. Results in Figure 2 (b), in turn, are as in Figure 2 (a) except that the inversion was done using only the samples #2, #3, and #4. This represents a more difficult and realistic case, since fewer mixtures are available.

When comparing Figure 2 (a) with Figure 2 (b), one can conclude that the method performs slightly better when more samples are used to estimate the true expression profiles – a result that was expected. Overall performance, however, is good in both cases. The estimated expression values for the pure colon cancer (RKO) are close to the mixture #1, as it should be since the mixture #1 corresponds to a measurement of the pure colon cancer. Similarly, the estimated expression values of the pure lymphocytes are close to the mixture #5 as well as to all of the reference samples (note that samples used in the reference channel (Cy5) are the same lymphocytes as the ones used in the mixtures). In Figure 1 and 2 (a), the most significant PCA component and the two most significant PCA components explain about 70.0% and 81.9% of the total variation in the data, respectively. For the reduced data, for which the results are shown in Figure 2 (b), a slightly smaller fraction of the variance is explained, namely about 67.3% and 81.2%. The results obviously depend on the optimality criterion for which we used the standard least squares. Less outlier sensitive results can be obtained with robust regression methods, such as the Huber estimator with the iteratively reweighted least squares implementation [12, 13] or median based regression methods [12, 14]. The robust methods provided similar global results, but improved results for some individual genes that contained one or more outliers.

Optimization of mixing percentages

In practice, the true mixing percentages are not known but must be measured by some means. Therefore, they are also likely to contain some error. So, in addition to inverting the mixing effect, it is also useful to simultaneously estimate the most likely value of the mixing percentages. This problem can be formulated as type of optimization problem, the details of which are shown in the Materials and methods Section. The proposed optimization scheme was applied to the heterogeneous microarray data. Since the heterogeneous samples #1 and #5 correspond to the cases where only colon cancer cells and lymph node cells are used, respectively, we may assume that α₁ = 1 and α₅ = 0. Thus, we only estimate the value of the three remaining mixing parameters. However, practically the same results are obtained when all the five mixing parameters are estimated. We found that the convergence of the above method is practically independent of the initialization in step 1. The convergence of the optimization method is illustrated in Figure 3 by showing the evolution of the value of the objective function. Parameters in Â⁽¹⁾ are initialized using the measured values shown in Table 3.

The obtained expression estimates for the pure colon cancer and lymph node samples, when all five heterogeneous samples are used in estimation, are shown in Figure 4. Again, the two most significant PCA components of all the heterogeneous samples, reference samples, and the estimated expression profiles of the pure colon cancer cells and lymphocytes are shown. It is instructive to compare these results with the ones shown in Figure 1. Because the heterogeneous samples are again located almost on a straight line, we use 1-dimensional visualization for the results. Figure 5 shows the obtained expression estimates in 1-dimensional PCA space.

Again, the estimated expression values for the pure colon cancer cells (RKO) are close to those from mixture #1, as it should be, since mixture # 1 corresponds to a measurement of the pure colon cancer cells. Similarly, the estimated values from the lymph node sample are close to those from mixture #5 as well as to all of the reference samples. In Figures 4 and 5 (a), the first PCA component and the first two PCA components explain about 69.2% and 81.7% of the total variation in the data. For the reduced data, for which the results are shown in Figure 5 (b), the fractions of variance explained are about 65.1% and 80.4%. Although the fraction of variance explained is slightly smaller than without the optimization of the mixing parameters, the optimized mixing parameters provide a better fit to the data.

Confidence intervals

Above we were only interested in estimating the expression values of the pure cell types. Often it is also useful to assess the confidence intervals of the obtained expression estimates. For that purpose, we consider two methods: one based on Gaussian approximation and the other using bootstrap. (For more details, see Materials and methods Section.)

For illustration purposes, Figure 6 shows the 90% estimated confidence intervals for a set of genes by pooling each of them with the 50 closest genes. The horizontal and vertical axes correspond to the fraction of lymph node cells and the normalized expression value, respectively. In other words, the different heterogeneous mixtures are placed on the x-axis according to the corresponding mixing fractions. The vertical lines at x = 0 and x = 1 expand over the maximum of the two confidence intervals. In most of the cases the two confidence intervals are in good agreement. The confidence intervals can be tightened by measuring more heterogeneous mixtures.

The proposed inversion methods for the sample heterogeneity were also tested on standard non-replicated microarray data by treating the replicated measurements for each gene as individual "genes." The obtained results were qualitatively similar with the ones shown above and only slightly more variable. In a similar fashion, we examined the effect of low quality replicates on the heterogeneity inversion. Slightly less variable results were obtained with a method [15] that detects and removes unreliable replicates prior to the averaging. A drawback of such unreliable spot detection is that, without any missing value estimation method, some of the genes will be excluded from further analysis.

Selection of the number of cell types

It is known that the heterogeneous mixtures used in our experiments consist of only two cell types. However, in general case, heterogeneous mixtures may contain an unknown number of cell types. In those cases, it is useful to assess the validity of the model (i.e., the number of cell types) as well. As introduced in Materials and methods Section, the linear mixing model can be extended to incorporate more than just two cell types. We use a general purpose cross-validation for model selection. In particular, we apply the so-called leave-one-out cross-validation (LOOCV) and test the one, two, and three cell type models. (For more computational details, see Materials and methods Section.)

Application of 'in silico microdissection' to detection of differentially expressed genes

In order to illustrate the above 'in silico microdissection' in practice, consider the following (hypothetical) experimental setting. Given the three middle mixture measurements (#2, #3, and #4), a goal is to identify a set of genes which are differentially expressed between the colon cancer and the lymph node samples. In a simple approach, often used in practice, the most heterogeneous sample would be discarded since it is measured to contain about 56% (resp., 44%) of colon cancer cells (resp., lymphocytes), thus giving no direct discriminative information about the underlying two samples. For illustration purposes, let us measure the expression difference of a given gene between these two samples using the fold-change, i.e., the expression value of the i th gene in the colon cancer sample, , is regarded as being differentially over-expressed (resp., under-expressed) if the ratio of to the expression value of the same gene in the lymph node sample, , is at least 2 (resp., smaller than 1/2). Of course, in practice, more sophisticated methods for detecting differential expression, including correction for multiple testing, should be used. However, for illustrative purposes, this example will suffice. Since only the heterogeneous samples are available, without any inversion of the mixing effect, one must compare the mixture measurement and . Figure 7 shows some example genes whose expression difference (i.e., the fold-change) between the two heterogeneous samples is within the given threshold (above 1/2 and below 2), but after the 'in silico microdissection,' the expression difference exceeds even a more stringent criterion (approximately 4-fold-change). The measured mixing percentages are used in the estimation (see Table 3). It is clear from this example that the proposed method is able to correctly detect differential expression even from heterogeneous samples, especially when the direct use of such samples may fail to find differential expression. Indeed, the conclusions we can draw based on the red stars are consistent with those that are based on the true homogeneous samples represented by blue squares in Figure 7.

As is evident from the example above, heterogeneity in the biological sample preparation can hinder further statistical analysis steps. Not only can the heterogeneity blur the identification of differentially expressed genes, it can also cause contrary effects. Presence of a considerable percentage of additional cell types can result in the identification of differentially expressed genes that may be unrelated to the biological question being studied. Similarly, irrelevant gene combinations can be discovered in the case of gene expression based classification. For an illustration, see [16] where the authors analyzed a colon cancer data set contaminated with muscle cells.

Although the microarray technology has been improved during the recent years, the measurements are still moderately noisy. The easiest and the most widely used approach for improving the measurement quality is to capture replicated measurements. This may become costly because each additional measurement requires an extra spot on the array, or an extra array. An alternative approach based on so-called composite microarrays was introduced in [17], where several different oligos representing different genes are printed on each spot. The multiplexing results in a mixing effect similar to the one introduced in this manuscript, and the phenomenon can be inverted to get the reconstructed expression values for single genes. The benefit is to obtain more replicated measurements without proportionately increasing the number of printed spots. Closely related ideas have also been introduced from an error-correcting microarray design point of view in [18]. The standard non-repeated microarray method does not tolerate "drop-outs": if a spot is badly corrupted and its intensity cannot be read, the expression value of the corresponding gene will be missed. Khan et al. showed that a certain amount of "drop-outs" can be recovered from the multiplexed samples, thus providing more error-resilient measurements. Following the methods developed in [17, 18], instead of multiplexing individual genes on spots, one may wish to multiplex different samples on arrays, thus allowing a fault-tolerant recovery of expression values in the case of corrupted array(s). As a future extension, one can also consider multiplexing both the genes on spots and the samples on arrays. Similar methods for inverting the sample heterogeneity have also been studied in the context of time-series gene expression measurements in [19, 20], where the fundamental mixing effect is not due to the different tissue types present in the sample, but due to the loss of synchrony of the cell population. It would be worthwhile to simultaneously study the sample heterogeneity and the loss of synchrony in the future.

Microarray production

RNA isolation, microarray production, and microarray hybridization were carried out as described previously in [21]. RNA from normal human lymph node was purchased from a commercial source (Stratagene, La Jolla, CA). Five μ g aliquots of total RNA from normal lymph node and RKO colon cancer cell line were reverse transcribed using Superscript II RT (Invitrogen, Carlsbad, CA) in conjunction with oligodT-T7 primers according to the manufacturer's suggested protocol. The second strand was synthesized using 10U E. coli DNA ligase (vendor), 40U E. coli DNA polymerase I (vendor), and 2U E. coli Rnase H (vendor). This reaction was stopped with EDTA and then cleaned with Qiagen's PCR Purification kit (Qiagen, Valencia, CA). The double stranded cDNA was then amplified by an in vitro transcription reaction (Ambion, Austin, TX) and cleaned with Qiagen's Rneasy kit (Qiagen, Valencia, CA). Each amplified cRNA sample was then quantitated using a Beckman DU640 spectrophotometer (Beckman, Fullerton, CA). Five μ g amplified cRNA from Stratagene's normal lymph node was labeled with Cy5 for each microarray hybridization. Mixtures of appropriate volumes of cRNA from normal lymph node and RKO were labeled with Cy3 in a reverse transcription reaction using Superscript II RT (see Table 3). Labeled samples were co-hybridized overnight at 60°C in a humidified incubator on a cDNA microarray containing 4704 human genes in duplicate produced in-house. The 4704 genes represent most of the known genes in the cDNA library we used to generate the microarrays. For the purpose of this study, the identity of the genes is not very important since we only study the general effect of sample heterogeneity. As the mixing effect is the same for all the genes, we expect to have similar results when the whole genome arrays are used. Slides were scanned with an LS-IV laser scanner (Genomic Solutions, Ann Arbor, MI). In total, five different heterogeneous mixtures were measured. The measured mixing percentages are shown in Table 3.

Preprocessing

The microarray data consists of five different heterogeneous mixtures of lymph node and colon cancer samples which are hereafter abbreviated as normal and RKO, respectively. (For more details, see Microarray production Section above and Table 3.) The gene expression data set was preprocessed as follows. The replicated background-subtracted signal intensities were averaged and log₂-transformed, and the dye-bias effect was corrected in the log₂-domain using the standard lowess smoothing-based normalization (see e.g. [22]) with smoothing parameter f = 0.7. Because the averaging effect (source of heterogeneity) takes place on the molecular level, the phenomenon must be modeled using the absolute expression values. Therefore, after the correction of the dye-bias, the data were transformed back to the original domain using the inverse of the log₂-transformation. Correspondingly, single-channel data were used for further analysis. In order to mitigate the between array variability, the data were further standardized for each array and the two channels separately.

Modeling sample heterogeneity

The two samples, RKO colon cancer cells and normal lymphocytes, are mixed at the extracted RNA level. Therefore, without any further verification, the model can be assumed to be linear. Lymphocytes were used because tumor tissues often contain infiltrating lymphocytes, especially in tumor metastases in the lymph nodes. Let

and denote the expression level of the i th gene in the colon cancer (RKO) and in the lymph node (normal) samples, respectively. Assuming only two different cell types are mixed, the sample heterogeneity is modeled by a simple linear model

where

denotes the expression value of the i th gene in the k th heterogeneous sample, and 0 ≤ α _k ≤ 1 denotes the fraction of the colon cancer cells in the kth mixture. It is worth noting that we use the same mathematical model for the sample heterogeneity as in [9–11]. Also note that in Equation (1) it is assumed that the expression level in RKO ( ) and normal ( ) is "fixed" and does not change between heterogeneous measurements. In other words, the measurements come from the same heterogeneous sample with different mixing fractions. In order to allow variation in the expression values between different samples/treatments/time points, the same model can be applied separately to each set of measurements from the other samples/treatments/time points. The same model can also be extended to more than two cell types (for more details, see Selection of the number of cell types Section below).

Inversion of sample heterogeneity

The first objective is to invert the mixing effect shown in Equation (1), that is, to obtain estimates for the expression values of the pure colon cancer cells and the pure lymphocytes. In practice, however, the measured expression values, y _i , include one or more sources of noise. By making some distributional assumptions, one could use standard model-based estimation methods. However, in order to avoid making additional modeling assumptions, we prefer to use a general purpose least squares method to estimate the gene expression levels corresponding to the pure samples.

Let the number of genes be n and assume that one has measured the expression values for K different heterogeneous mixtures. Thus, one has measurements , 1 ≤ i ≤ n, 1 ≤ k ≤ K. Let us also assume for now that the mixing percentages are known or have been measured. For the i th gene the sample heterogeneity can be expressed as (excluding all noise terms)

When including all n genes, the above model can be rewritten as

where 0 denotes the K-by-2 zero matrix. Let the block matrix in Equation (2) above be denoted as à Assuming the column rank of A is full, the well-known least squares solution is given by

where

. Due to the structure of the matrix Ã, the least squares solution can be obtained gene-wise as .

The Gauss-Markov theorem says that the standard least squares solution is indeed the best linear unbiased estimate if the noise in the measurements is additive and i.i.d. with constant variance. However, a common observation is that the homoscedasticity does not always hold for microarray data, but instead, the noise variance depends on the underlying signal intensity [23, 24]. Such heteroscedasticity may decrease the power of the inversion method shown in Equation (3). Fortunately, the structure of the matrix à ensures that the inversion can also be performed for each gene separately. Consequently, it is not necessary for the homoscedasticity to hold globally. Indeed, all we need to assume is that the noise variance is approximately constant for each gene separately.

Also note that, in this two cell type model, no prior knowledge about the expression values of either of the two cell types is needed since the method estimates the expression values for both of the two cell types. The same is also true for more general models including more cell types, assuming the model is sufficiently over-determined (see also Selection of the number of cell types Section below).

Optimization of mixing percentages

In practice, the mixing percentages must be measured by some means. Therefore, they are also likely to contain some error. So, overall, one would like to estimate not only the expression values for the pure cell types but also the most likely value of the mixing percentages.

Assuming the model in Equation (2) is sufficiently over-determined, the mixing parameters can be adjusted computationally, too. Let us again consider the case where only two different cell types are mixed. Note that K denotes the number of different heterogeneous mixtures measured. Therefore, the regression matrix à in Equation (2) has only K free parameters. Since the number of expression values to be estimated is 2n, the total number of free parameters in Equation (2) is 2n + K. The number of equations in Equation (2) is Kn. Hence, the model is over-determined if Kn > 2n + K, which, for a fixed n ≥ 3, holds if K > 2. (Note that in our case we have measured five different heterogeneous mixtures, i.e., K = 5.) As above, no assumptions on the noise distributions are being made and we use the least squares method. This results in the following optimization problem

A similar optimization problem was also introduced in [11].

Because the objective function in Equation (4) above is minimized over both à and x, the objective function is not linear in the parameters anymore and, therefore, cannot be solved as in Equation (3). In general, any iterative optimization method can be used to get a solution. Iterative methods usually become inefficient/unstable as the number of parameters to be optimized increases. In this case, the number of free parameters in à and x is 2n + K. Therefore, we use a two-step approach in the optimization. In the first step, given proper initial value for Ã, the least squares solution for x is found using Equation (3). In the second step, the mixing percentages are optimized in the least squares sense (subject to the constraints 0 ≤ α _k ≤ 1 for all 1 ≤ k ≤ K) using the previously found value for x. These two steps are then repeated, essentially resulting in a type of expectation-maximization (EM) approach. A similar iterative procedure was also proposed in [11], except with different constraints. Note that when Equation (4) is minimized over Ã, given the value of x, the optimization problem is again linear in its parameters. Assuming the constraints are not violated, the standard equation (similar to the one in Equation (3)) can be applied. If that is not the case, then any general-purpose constrained optimization method may be applied. Let (resp. Â^(j)) denote the value of x (resp. Ã) after the j th iteration. Details of the algorithm are shown in Figure 8. Clearly, at each iteration of steps 2 and 3, the value of the objective function is decreased. Thus, a minimum will be found.

Confidence intervals

It is useful to assess the confidence intervals of the obtained expression estimates. As explained above, the Gauss-Markov theorem is applied gene-wise that greatly alleviates the issue of heteroscedasticity. Should the noise variance σ² be constant, then the variance of the estimated expression values would be . Due to the special structure of the matrix à (i.e., the gene-wise inversion of the mixing effect), the variance of the estimated expression values for the i th gene can be expressed as

where

is the noise variance for the i th gene. A straightforward way of obtaining an estimate of the variance is to compute the sample noise variance for each gene and then apply Equation (5) to get . That would result in somewhat sensitive variance estimates since there are only K = 5 error residuals associated with each gene. A better alternative is to pool genes which have approximately the same average expression value and then compute the sample noise variance from the error residuals of the pooled genes. Although we do not assume Gaussian noise distribution, we can resort to the Gaussian approximation when computing the confidence intervals. For example, using the Gaussian approximation, the 1 - 2α confidence interval for estimated expression value of the i th gene in the colon cancer cells is , where Φ^-1(·) is the inverse of the standard normal cumulative distribution function and denotes the (1,1) element of the estimated variance matrix (similarly for the lymph node sample: ). Alternatively, the confidence intervals can be obtained using the non-parametric bootstrap framework [25]. Here we consider the method in which one re-samples the error residuals with replacement (within the set of pooled genes) and computes the confidence intervals directly from the α and 1 - α percentiles of the bootstrap distribution of the expression estimates.

Selection of the number of cell types

Although it is known that only two cell types are mixed in our experiments there may be other experimental settings where the number of cell types may be unknown. Then it is useful to assess the validity of the model as well. As was mentioned above, the linear mixing model can be extended to incorporate more than just two cell types using a straightforward extension:

, where denotes the expression value of the i th gene in the j th cell type, and 0 ≤ ≤ 1 denotes the fraction of the j th cell type in the k th mixture. The mixing percentages must also satisfy for all k. The significance of different 'regression coefficients' could be tested using standard regression-based statistical tests. Since those tests apply only to Gaussian noise we recommend using a general purpose cross-validation for model selection (see e.g. [26]). Here we consider the leave-one-out cross-validation (LOOCV) and test the one, two, and three cell type models. Thus, each heterogeneous sample is left out from the training data at a time, the regression coefficients are estimated based on the remaining four samples, and the model is then tested on the sample which was left out from the training data set.