Robust computational reconstitution – a new method for the comparative analysis of gene expression in tissues and isolated cell fractions

Background Biological tissues consist of various cell types that differentially contribute to physiological and pathophysiological processes. Determining and analyzing cell type-specific gene expression under diverse conditions is therefore a central aim of biomedical research. The present study compares gene expression profiles in whole tissues and isolated cell fractions purified from these tissues in patients with rheumatoid arthritis and osteoarthritis. Results The expression profiles of the whole tissues were compared to computationally reconstituted expression profiles that combine the expression profiles of the isolated cell fractions (macrophages, fibroblasts, and non-adherent cells) according to their relative mRNA proportions in the tissue. The mRNA proportions were determined by trimmed robust regression using only the most robustly-expressed genes (1/3 to 1/2 of all measured genes), i.e. those showing the most similar expression in tissue and isolated cell fractions. The relative mRNA proportions were determined using several different chip evaluation methods, among which the MAS 5.0 signal algorithm appeared to be most robust. The computed mRNA proportions agreed well with the cell proportions determined by immunohistochemistry except for a minor number of outliers. Genes that were either regulated (i.e. differentially-expressed in tissue and isolated cell fractions) or robustly-expressed in all patients were identified using different test statistics. Conclusion Robust Computational Reconstitution uses an intermediate number of robustly-expressed genes to estimate the relative mRNA proportions. This avoids both the exclusive dependence on the robust expression of individual, highly cell type-specific marker genes and the bias towards an equal distribution upon inclusion of all genes for computation.

(k) s ]. In addition, the standard deviations of the relative proportions are calculated from the inverse of the second derivatives of E[D (k) s ] [1]. The good correspondence to the numerical results for the patient data shown in Figure 2 of the main text corroborates the suitability of the proposed method for the determination of the relative mRNA proportions. The third part of this section includes graphs of the empirical objective function D (k) s (trimmed sum) illustrating its dependence on the number of included genes k.

A.1 Convergence of the relative proportions
The model used for showing the convergence of the computed relative mRNA proportions is strongly simplified in that the difference δ s = log S −log S * (with S = pMM +pFF +pNN and S * = p M M +p F F +p N N ) is approximated by in whichm = logM , m = log M , etc. This approximation was necessary because the full form of δ s was not amenable to an analytical treatment. It can be justified in a statistical sense in that the measured log-transformed expression profiles s, m, f and n are all approximately normally distributed, which is also true for a linear combination of normally distributed random variables. The problem was further simplified in that all auto-correlation coefficients (e.g., rm m ) were set equal to r, all cross-correlation coefficients (e.g., r mf ) were neglected, and all standard deviations were assumed to be equal to σ, because only the principal mechanisms were of interest. Under these assumptions, the difference δ s is normally distributed with the standard deviation [2,3]

A.2 Model for Robust Computational Reconstitution
A realistic model for Robust Computational Reconstitution must allow for the fact that the correlation coefficient r is generally different for each gene. This was done by multiplying the Gaussian probability density function (PDF) of δ s having standard deviation according to equation (A.2) (and thus a specific value for r) by an empirically determined piecewise linear PDF for the correlation coefficient r (Figure A.1) and subsequent marginalization with respect to r. The resulting modified PDF for δ s was then used for calculating the expectation values E[D (k) s ] of the trimmed sum D (k) s in Equation (2) of the main text. This was done by summing up the expectation values of the lowest k values of the absolute difference ρ = |δ s | in a random experiment with N (number of genes) trials. The expectation values for these k lowest values were calculated using order statistics [4,5]. Let F (ρ) be the cumulative distribution function (CDF) of ρ. Then, the CDF of the ith smallest value ρ i (ith order statistic) is given by Using one partial integration, the expectation value of ρ i was calculated according to with f i (ρ i ) denoting the PDF of ρ i andρ being an appropriate cut-off value, which was determined according to 1 − F (ρ) N < 10 −8 . The expectation value of D (k) s = k i=1 ρ i was then calculated by summing up the first k expectation values E[ρ i ]. Finally, the relative proportions p M , p F , and p N that minimize E[D (k) s ] were calculated as a function of the number of included genes k. The results shown in Figure A.2 are similar to the numerical results obtained for the experimental data of patient 2 ( Figure 2 of the main text) but give the best estimates for the relative proportions at the lowest number of included genes. This is because the expectation value E[D (k) s ] is a smooth function of the proportions and the minimum can accurately be determined. For a specific data set, however, a multitude of equivalent local minima is scattered across the whole domain (triangle 0 ≤ p M + p F ≤ 1) implying an equal distribution of the relative proportions. Small changes in the PDF of the correlation coefficient r near r = 1 have a strong influence on the maximum attainable precision of the estimated relative proportions. This is demonstrated by the differences between the dashed and the solid curves in Figures A.1 Figure 2 of the main text, however, they do not converge towards an equal distribution for a small number of included genes because the expectation value E[D (k) s ] is a smooth function of p M and p F (see text). The standard deviations drop earlier for the same reason. Solid line: medium density assumed at r = 1, dashed line: zero density assumed at r = 1 (see Figure A.1). The number of included genes ranges between 25 and 1253 (one 10th of the 12533 genes on the chips) because summation and integration according to Equations (A.4) and (A.5), respectively, was done numerically.

A.3 Empirical objective function of trimmed regression
In the present study the trimmed sum D (k) s in Equation (2) of the main text is used as the objective function for determining the relative mRNA proportions of the tissue samples of each patient. The figures shown in this section illustrate the dependence of D (k) s on the number of included genes k, which is essential to the method of robust computational reconstitution. The objective function turned out to be very similar for L 1 (least modulus) and L 2 (least squares) regression. Generally, L 1 regression was used in the present study because it was expected to be more robust with respect to y-outliers. For comparison, the results for patient 2 using L 2 regression are shown in Figure A

B Experimental Materials and Methods
Pieces of synovial tissue obtained from 3 patients fulfilling the American Rheumatism Association criteria for rheumatoid arthritis (RA, patients 1-3) [6] and 3 patients with osteoarthritis (OA, patients 4-6) [7], were placed in RNAlater (Ambion Inc., Austin, TX, U.S.A.) for 24 h at room temperature and then frozen immediately at -70 • C. The different synovial cell populations were prepared as previously described [8]. Briefly, the tissue samples were minced, digested with trypsin and collagenase P, and the resulting single cell suspension cultured for 7 days. During this primary culture, non-adherent cells were removed by medium exchange on day 1. After 7 days, synovial fibroblasts were purified (98% purity) by removal of macrophages using Dynabeads R M-450 CD14 and the Dynal magnetic particle concentrator R . Non-adherent cells, fibroblasts and Dynabeads-coupled macrophages (>95% purity) were immediately lysed in RLT buffer (Qiagen, Hilden, Germany) and frozen at -70 • C. Total RNA was isolated using the RNeasy Kit (Qiagen) according to the supplier's recommendation. RNA probes were labeled according to the supplier's instructions (Affymetrix, Santa Clara, CA, USA). Analysis of gene expression was carried out using HG-U95A (patient 1) and HG-U95Av2 (patient 2-6) chips. Hybridization and washing of the gene chips was performed according to the supplier's instructions (http://www.affymetrix.com). Microarrays were analyzed by laser scanning (Hewlett-Packard Gene Scanner). The synovial tissue chip of patient 1 showed about twice as much noise as the other chips according to the operating figures provided by Affymetrix' MAS 5.0 software. Therefore, the hybridization was repeated twice (using HG-U95Av2 chips), the second repeat having been more stringently washed. Finally, the first repeat was selected for analysis because its operating figures were most compatible with the remaining chips. The purity of the isolated cell fractions was confirmed by a strong down-regulation of specific marker genes for macrophages, fibroblasts, and T-cells (a component of non-adherent cells) in the respective other cell fractions. Two own mixing experiments were performed, in which mRNA preparations of the isolated cell fractions of patient 2 were mixed according to two sets of relative mRNA proportions (that are representative for the range of computed values) and were then hybridized to a HG-U95Av2 chip. Immunohistochemical analysis was performed on cryostat sections of RA synovial membranes (fixed with acetone for 10 min at 4 • C) using the following monoclonal antibodies (mAbs): AS02 (CD90; Dianova GmbH, Hamburg, Germany) and 3-2B12 (prolyl-4-hydroxylase; Merck Biosciences, Schwalbach, Germany; both fibroblast markers); CLB CD14B (CD14; Janssen Biochimica, Beerse, Belgium; monocyte/macrophage marker); UCHT-1 (CD3; Dako, Hamburg, Germany; T-lymphocyte marker); HIB19 (CD19; BD, Heidelberg, Germany; B-lymphocyte marker); AT13/5 (CD38; Dako; plasma cell marker); 4F9 (van Willebrand Factor; Beckman Coulter, Krefeld, Germany; endothelial cell marker). The primary mAbs, diluted in phosphate-buffered saline (PBS)/1% horse serum, were added to the tissue sections for 30 min in a humid chamber at room temperature following pre-incubation with 5% normal horse serum in PBS for 20 min. A peroxidase-coupled rabbit anti-mouse antibody (Dako; 1:30 diluted in PBS/1% horse serum) was then added for 30 min. The peroxidase was revealed under visual control for approx. 5 min using diaminobenzidine (0.5 mg diaminobenzidine in 1 ml PBS containing 30 µl H 2 O 2 ). The slides were washed, counterstained with haematoxylin and covered with Aquatex (Merck). Isotype-matched control mouse mAbs, used at identical concentrations as the primary antibodies, gave no positive results.
For the assessment of the cellular composition, one to three different synovial tissue sections from each patient were analyzed. The cellular composition was separately assessed in the following compartments: lining layer, diffuse infiltrates, lymphoid aggregates, and stroma. The area of the individual compartments was measured in each section using the software AxioVision 4.2 (Carl Zeiss Vision GmbH, Hallbergmoos, Germany) and the relative percentage of cells showing a positive reaction for each marker in each compartment was then visually estimated in 10% steps. In addition, the cellular density of the respective compartment was determined by counting the cell nuclei in two to five 150 µm 2 areas. The percentage of cells positive for each marker in each compartment was then corrected for: i) the relative area of the compartment in the tissue section; and ii) the cell density of the compartment.

C Comparison of MAS-S and MAS-C probe set summaries
The MAS-S probe set summaries were generally preferred to the MAS-C, RMA, and MBEI summaries because they showed a better agreement between the mRNA proportions computed with and without additional stepwise local regression normalization. This is demonstrated by comparing the results obtained from MAS-S and MAS-C for patients 4 and 1 in Figures C.1 and C.2, respectively. The computed relative mRNA proportions are displayed in Table C summaries normalized by a symmetric trimmed mean (t) and quantile normalization (q). The reconstituted tissue profile is either not further normalized (dashed) or normalized by stepwise local regression (lr ) (solid). Similar results were obtained for the other three normalization methods (t, t clr, t c; for notation see Table 1 of the main text).   (1)

D Experimental chip quality
The differences in the computed mRNA proportions resulting from different chip evaluation methods, i.e. probe set summaries and chip normalizations, could be caused by varying experimental chip quality. In order to address this hypothesis, the MAD of the computed mRNA proportions (as defined in Table 2 of the main text) was correlated with several relevant chip operating figures provided by the Affymetrix 'Expression Report'. The operating figures estimate the background and noise level, as well as the degree of mRNA degradation for each chip. After scaling by the 'Scale Factor (SF)' of the respective chips, the operating figures were summarized for each patient by taking a weighted average across the four chips per patient. With a few exceptions, the MAD of the calculated mRNA proportions was positively correlated with the noise and background level ( Figure D.1). The correlation coefficients ranged from 0.60 to 0.93. However, due to the small number of patients the results strongly depended on individual patients (inhomogeneity correlation), in particular on patient 1 (original tissue chip). The methodological scatter and the correlation with the operating figures were considerably reduced when only the four different normalization methods for the MAS-S summaries were considered ( Figure D.2). The correlations became negative (correlation coefficients down to -0.89), when the additional stepwise local regression normalization was applied to the differently normalized MAS-S summaries, indicating that local regression can make the computed proportions more homogeneous in the presence of noise.
The correlation between the MAD of the computed relative proportions and mRNA degradation (as quantified by the 3'/5' expression ratios of GAPDH and β-actin) was generally weak, quite unsystematic and mostly negative. This suggests that not mRNA degradation, but rather experimental noise and background are the main parameters affecting the determinability of the tissue mRNA proportions.

F Cell type-specific marker genes
Stable cell type-specific marker genes allow to determine the mRNA proportions from their expression value in the synovial tissue and the respective isolated cell fraction (e.g., p F = s i /f i for fibroblasts and marker gene i). In the present study, the correspondence between the relative mRNA proportions calculated from Robust Computational Reconstitution and from marker genes were reasonably consistent only for fibroblasts (

G Regulated and robustly-expressed genes
Regulated and robustly-expressed genes were selected according to large and small differences, respectively, between their gene expression in the tissue and the respective isolated cell fractions (reconstituted tissue profiles) for all patients. Six different statistical methods (described in Section 2.5 of the main text) were used for this purpose. Each individual method was independently applied for gene ranking. A number of top-ranking genes were defined according to thresholding values (Table G.1) and subsequently screened for pathophysiological relevance. The genes identified as regulated and robustly-expressed, along with their particular identification methods, were finally listed in Tables G.2     The average overlap is 79.2 for regulated and 76.7 for robustly-expressed genes. The largest overlap for regulated genes is between µ-test and MAD-test, that for robustly-expressed genes is between SAM and t-test2.  SAM, t-test1, and t-test2 have a tendency to select highly-expressed regulated genes (low p-value range) and weakly-expressed robustly-expressed genes (high p-value range). The selection by V&S and µ-test appears to be almost unbiased with respect to M. The MAD-test has a clear preference for highly-expressed genes at the upper p-value end (robustly-expressed genes). This may reflect the smaller differences A = log 10 S − log 10 S * observed for highly-expressed genes (large M) in socalled M-A-plots (see also Figure H.3). Moving averages and standard deviations were calculated using a symmetric 200-genes window.  , with r = mean (log (S/S * )), S measured tissue gene expression, and S * reconstituted tissue gene expression. Mean expression < S > = exp (mean (log S)), < S * > = exp (mean (log S * )). P: number of present calls for the 6 synovial tissue chips, P * : number of present calls for the 18 isolated fraction chips. Data: trimmed mean normalized MAS-S probe set summaries. Number of patients: n = 6. Gene ranking was exclusively performed using the statistical methods. The values of < fold >, < S >, < S * >, P , and P * are only given as additional information.
gene name symbol probe set methods < fold > < S > < S * > P P *   Table G

H Probabilistic model for cell type-specific gene expression
The gene expression values of each individual cell type in the tissue can be estimated from the gene expression values of the isolated cell fractions and the whole tissue based on a suitable, fully specified probabilistic model, i.e. a model in which all distributional parameters are settled. In this section a multivariate normal distribution model is assumed for this purpose. Clearly, the cell type-specific mean expression values in the tissue as well as the respective correlation coefficients between tissue and isolated cell fractions cannot be estimated directly because the cell typespecific gene expression in the tissue is not available. However, these distributional parameters (mean and correlation coefficients) could possibly be approximated as an universal function of expression strength. This general approach of substituting statistical parameters using functions of expression strength can be tested for the standard deviations because they can be calculated in two alternative ways: either data-based using the expression values measured for the isolated cell fractions or by substitution, using a function of expression strength. Although the results of both methods were largely consistent, a number of contradictory results were observed (Table H.1). Consequently, the present substitution approach cannot fully compensate for the lack of experimental microdissection data and thus the cell type-specific gene expression values in the tissue cannot reliably be estimated. The probabilistic model assumes the log-transformed expression of each gene in each cell type to be a normally distributed random variable (subsequently denoted by lower case letters, omitting the gene index i introduced in Equation (3) with µ denoting the vector of mean values and Σ being the covariance matrix: Σ αβ = r αβ σ α σ β with r αβ being the correlation coefficient between α and β, and σ α and σ β representing the respective standard deviations. In order to ensure the computed tissue expressions(x) = log p M expm + p F expf + p N expn to be close to the measured value s, the PDF in (H.1) was regularized by multiplication with the Gaussian Analytical expressions for the expectation values and standard deviations of the cell type-specific changes in gene expression could be obtained by linearizings(x) with respect to the differences δ α =ᾱ − α for α = m, f, n (i.e. substitutingm = m + δ m and expanding with respect to δ m , etc.). Assuming, for notational convenience, equal mean values in tissue and isolated cell fractions, equality of means and expression values among the isolated cell fractions, and considering the cell type-specific correlation coefficients only (i.e. r αβ = 0 except for α = β and α =β), the expectation values E[δ α ] and standard deviations S[δ α ] of the expression differences read in which δ s = s − s * denotes the observed difference in gene expression, p · w = α p α w α is the scalar product between p and w, and p m = p M (M/S * ), etc. The conditions(x) = s is exactly met forσ s = 0 implying δ s = α p α E[δ α ]. Clearly, the numerator of E[δ α ] equally depends on the fraction p α , one minus the squared correlation coefficient r 2 αα , and the variance σ 2 α , i.e. the correlation coefficients rᾱ α need to be known or estimated.
For a better identification of simultaneous changes in groups of patients the above approach was generalized to multiple patients having identical changes δ α in gene expression. Assuming the patients to be statistically independent the total PDF is the product of the regularized PDFs for each patient introduced in Equations (H.1) and (H.2) In Equation (H.4) the cell type-specific gene expression values in the tissue were framed asᾱ p = α p +δ α for each patient and α = m, f, n . Also, Ω P summarizes the distributional parameters for all P patients. The standard deviationŝ σ sp of the regularizing Gaussians were obtained by approximating σ s = σ s (s) as a function of the mean expression level s using linear regression (Figure H.1) and multiplying the function values σ s (s p ) by a scaling factor, presently 0.1. Similarly, the cell type-specific correlation coefficients rᾱ αp were estimated by r ss * (α p ), i.e. by generalizing the dependence of the correlation coefficient r ss * between the measured and reconstituted tissue profiles (across patients) on the expression level to the individual cell types (Figure H.2). For assessing the reliability of the estimation method for rᾱ αp the cell type-specific standard deviations σ αp were determined in two alternative ways. On the one hand, they were directly calculated from the respective isolated cell fraction (across patients, thus being equal for all patients). On the other hand, they were determined using expression-level-based regression, analogous to the estimation of rᾱ αp described above (cf. Figure H.1).
Weakly-expressed genes show a tendency to be upregulated (and vice versa) (Figure H.3). This was incorporated into the model by assuming the differences between the mean values in the tissue and the respective isolated cell fractions to depend on the expression level (again via linear regression). However, the results obtained were very similar to the results reported in Table H.1 for equal means, indicating that this trend is too week to have a distinct effect on the results. Table H.1 shows the results for some arbitrarily selected genes. The cell type-specific assignment of the changes in tissue gene expression δ s sometimes was contrary between the two alternative ways of determining the standard deviations σ αp described above. Thus, in view of Figure H.2, which clearly shows the weak correlation between expression level and correlation coefficient, a reliable cell type-specific assignment of the changes in gene expression (tissue versus isolated cell fractions) appears to be impossible given the present data.      . Although the calculated cell type-specific fold changes agree in general, they can be reversed between the two methods, e.g. for the transforming growth factor, beta receptor III, and the bone morphogenetic protein 4, or can be markedly different, e.g. for interleukin 11 (numbers typed in bold face). The fold changes are based on the numerically calculated modus (most probable value) of the non-linearized PDF in Equation (H.4). The standard deviations of the fold changes (data not shown) depend on the choice of the scaling factor for the standard deviationσ sp of the regularizing Gaussians and were below 10% for all values shown (scaling factor 0.1). Expectation values and standard deviations were numerically calculated using a closed extended integration formula with the same order as Simpson's rule [9].