DNA methylation arrays as surrogate measures of cell mixture distribution

Statistical methods

Let Y_0h be an m × 1 vector of methylation assay values, e.g. average beta values from an Infinium bead-array product corresponding to a purified blood sample consisting of a homogenous cellular population (e.g. monocytes or granulocytes), with the qualitative characterization of cell type (among d₀ such types) indicated by a d₀ × 1 covariate vector w_h. Here, h∈{1,…,n₀}, where n₀ is the number of specimens and the m individual values correspond to CpG sites on a DNA methylation microarray, possibly pre-selected to correspond to putative DMRs for distinguishing different cellular types. Correspondingly, let Y_1i be an m × 1 vector of methylation assay values for the same CpG sites (in the same order) as Y_0h, but corresponding to a heterogeneous mixture of cells (e.g. peripheral whole blood) from a human subject. Here, i∈{1,…,n₁}, n₁ is the number of target specimens, and z_1i is a d₁×1 covariate vector representing phenotypes or exposures corresponding to the subject, e.g. d₁ = 2 for a simple case/control study without confounders. Our goal is to understand the associations between Y_1i and z_1i in terms of associations between Y_0h and w_0h, i.e. to infer changes in mixtures of cell types associated with phenotypes or exposures, using DNA methylation as a surrogate measure of cell mixture. Thus, we have two data sets, S₀ = {(Y₀₁,w₁),…,(Y_0n0,w_n0)}, the set of data from “purified” cell samples effectively representing external validation or gold-standard data, and S₁ = {(Y₁₁,z₁),…,(Y_1n1,z_n1)}, representing surrogate data collected from a target population. To this end, we posit the following linear models:

(1)Y0h=B0w0h+e0hY1i=B1z1i+e1i,

where B₀ and B₁ are, respectively, m × d₀ and m × d₁ matrices and e₀ and e₁ are error vectors. For simplicity we assume a one-way ANOVA parameterization for w, though in the Additional file1 we describe slight generalizations to account for design complications met in practice. We also assume a reasonable regression parameterization for z, including an intercept, and for convenience, denote the first column of B₀ as μ₁, the m × 1 intercept. The error vectors e₀ and e₁ may reflect independence among arrays h and i, or else may have more complex random effects structure accounting for technical effects or biological replication; however, their substructures are incidental to this analysis, with the exception of the fine details of the bootstrap procedure proposed below.

To implement a surrogacy relation, we propose the following linking regression model:

(2)B1=1mγ0T+B0Γ+U,

where Γ is a d₀ × d₁ matrix that summarizes associations between the rows of B_0j and B_1i and U is a matrix of errors. Substituting equation (2) into (1), writing B₀ = (b₀₁,…,b_0d0) explicitly in terms of its columns and writingΓT=(γ1,…,γd0), it follows that

(3)Y1i=∑l=0d0b0l(γlTz1i)+(1mγ0T+U)z1i+e1i.

To impart a biological interpretation, we assume that the DNA assayed in S₁ arises as a mixture of DNA from cell types profiled in S₀, with mixture coefficients whose population averages, conditional on z, are{ω1(z),…,ωd0(z)}, so that

(4)E(Y1i|z1i=z)=ξ(z)+∑l=1d0b0lωl(z),

where the m × 1 vector ξ^(z) represents cell types excluded from consideration among the purified samples in S₀, or else non-cell-specific methylation, including alterations at the molecular level in the maintanence of DNA methylation patterns themselves (possibly exposure related, age, or disease related). It follows from (3) and (4) that the mixture coefficients are recoverable from Γ,ωl(z)=γlTz1i, provided ξ^(z) is orthogonal to the column space of B₀. As we discuss in detail in the Additional file1, bias can arise if differences in ξ^(z) between distinct values of z have nonzero projection onto the column space of B₀, although the magnitude of anticipated biases can be assessed through sensitivity analysis.

It is possible to assign interpretations to the components of variation in (3). Let SS_o represents overall variability in Y_1i, i.e.SSo=∑i=1n1∥Y1i−μ¯1∥2, whereμ¯1=E(Y1i). From multivariate probability theory it is straightforward to show that SS_o = SS_e + SS_v + SS_u, whereSSe=∑i=1n1∥e1i∥2,SSv=∑i=1n1(z1i−z¯1)TΓTB0TB0Γ(z1i−z¯1), andSSu=∑i=1n1{(z1i−z¯1)TUTU(z1i−z¯1)+m(z1i−z¯1)Tγ0γ0T(z1i−z¯1)}. SS_e measures variation unexplained by the covariates z_1i, presumed to represent a combination of technical noise and unsystematic biological heterogeneity. SS_v measures variability explained by mixtures of profiles in the set S₀, while SS_u measures variability in systematic biological heterogeneity that nevertheless remains unexplained by mixtures of profiles in S₀, presumably due to some process other than differences in mixtures of cell types. Thus we propose two partial coefficient of determination measures:R1,02=SSv/SSo, which represents the proportion of total variation in S₁ explained by S₀, andR1,12=SSv/(SSo−SSe), which represents the proportion of systematic variation in S₁ explained by S₀. Note thatR1,12 is poorly defined when SS_o≈SS_e.

Estimation procedes by applying an appropriate linear model, e.g. ordinary least squares, linear mixed effects models[16], limma[17], or surrogate variable analysis[18,19], to obtain estimatesB^0 andB^1. Estimates of γ₀ and Γ are then obtained by projectingB^1 onto the column space ofB~0=(1m,B0), as described in detail in the Additional file1. Standard errors can be obtained in one of three ways. The simplest estimator, SE₀, is the “naive” estimator from simple least-squares theory, ignoring the fact thatB^0 andB^1 are estimates, i.e. potentially variable. To account for variation in estimatingB^1, a simple alternative is to use a nonparametric bootstrap procedure. For each bootstrap iteration t, we sample with replacement from S₁ (or sample errors in a manner consistent with a hierarchical experimental design) to obtainS1(t), producing bootstrap estimatesB^1(t) from which “single-bootstrap” standard errors SE₁ are computed. Finally, it is possible to account for variation in estimating _B0by also bootstrapping S₀; because of potentially small sample sizes n₀, we propose using a parametric bootstrap. A“double-bootstrap” standard error estimator, SE₂, is computed from these two sets of bootstraps. The double-bootstrap has the additional benefit over the single-bootstrap, in that it can be used to assess bias due to measurement error (variability) inB^0. Estimation details are provided in the Additional file1, as are the results of simulation studies.

Beyond bias due to measurement error, which is easily corrected using the double-bootstrap procedure, there are additional sources of potential bias. For example, consider a univariate z_1i representing case/control status, whereδ≡ξ(1)−ξ(0)=B0α for some d₀ × 1 vector α ≠ 0; i.e. δ is the mean difference in DNA methylation between a case and control, contributed by cell mixtures that remain uncharacterized or non-cell-specific methylation. In such a situation, there will be a bias equal to α in estimating the mixture differences. The Additional file1 provides a detailed analysis of such biases, and proposes a sensitivity analysis procedure for assessing the magnitude of possible bias in a given data set.

While the focus of this paper is analysis of population data, it is possible to use S₀ to predict distribution of leukocytes in a single sample having DNA methylation profile Y^∗. Equating the intercept term of B₁ in (1) with Y^∗ and applying (2), we obtain mixing proportion estimatesΓ∗=(B~0TB~0)−1B~0TY∗. Estimates can be further refined with the use of quadratic programming techniques[20], restricting the components of Γ^∗,γ_l^∗ ≥ 0, in minimizing∥Y∗−B~0Γ∗∥2 with respect to Γ^∗. Such individual projections of methylation profiles on the column space spanned by S₀ facilitate the application of the fundamental ideas proposed above to individual, clinically-based diagnostic procedures. Note, however, that DNA methylation arrays are typically focused on the comparison of methylated to unmethylated CpG dinucleotides, not quantifying actual amounts of DNA. Therefore, information on cell mixtures from DNA methylation is limited to distributions, not actual counts, as one might obtain from flow cytometry. Finally, we remark that it is possible to model z_1i directly as a function of mixture coefficients Γ^∗ obtained individually via the constraint γ_l^∗ ≥ 0, but the inferential implications are less clear, and we view the proposed approach for populations as more statistically robust.

Implementation

We describe several examples using existing methylation data sets as benchmarks for validating the proposed method, in order to demonstrate its clinical or epidemiological utility. First we describe the validation data set S₀ used in all examples. Next we describe a laboratory reconstruction experiment, which validates our fundamental proposition that DNA methylation retains substantial information about cell mixtures. Finally we describe the results of applying our methodology to several different target data sets S₁. For the head and neck cancer and ovarian cancer data sets, from which bead chip data were available, a linear mixed effects model with a random intercept for bead chip was used to estimate the corresponding row of B₁. For the remaining data sets, no bead chip data were available; consequently, ordinary least squares was used. 250 bootstrap iterations were used for each example and each of the two bootstrap methods of standard error estimation.

Validation data

All data analyses involve DNA methylation data obtained by the Infinium HumanMethylation27 Beadchip Microarrays from Illumina, Inc. (San Diego, CA). We used a subset of m = 100 CpG sites on the array, selected as described below. In all of our examples, S₀ consisted of 46 white blood cell samples, de-identified specimens that were not subject to human subjects review by an institutional review board (IRB). The sorted, normal, human, peripheral blood leukocyte subtypes were purchased from AllCells^Ⓡ, LLC (Emeryville, CA) and were isolated from whole blood using a combination of negative and positive selection with highly specific cell surface antibodies conjugated to magnetic beads; materials and protocols were obtained from Miltenyi Biotec, Inc. (Auburn, CA). These 46 samples are summarized in Table1 and depicted by the clustering heatmap in Figure1. Note that T lymphocytes that express CD4 or CD8 constitute over 95% of the T cell class, and that the pan-T cell type was further refined to CD4+, CD8+, and “other” Pan-T cells subtypes. In summary, the covariate vector w_hconsisted of indicators for five cell types and another two indicators for CD4+ and CD8+ T cell subtypes. A generalization of the one-way ANOVA parameterization assumed above for w_h, described in the Additional file1, was necessary to account for the ambiguous status of some Pan-T cells. For each CpG site, a linear mixed effects model with a random intercept for bead chip was used to estimate B₀; 27 additional whole blood control samples (replicates from the same individual) were used to assist in estimating chip effects, since otherwise the data set would have been sufficiently sparse to risk confounding between cell type and chip. These “array controls” were indicated with an additional term in w_0h. For each CpG site, a linear mixed effects model with a random intercept for bead chip was used to estimate the corresponding row of B₀and B₁. From S₀, F statistics (described in the Additional file1) were computed and used to order each of the 26,486 autosomal CpGs by decreasing level of informativeness with respect to blood cell types. As described in the Additional file1, we determined that maximum informativeness was provided by the top m = 100 − 300 CpG sites, with m > 300 reflecting diminishing returns from adding additional CpGs. Therefore, we chose a moderately low value in this range, m = 100, consistent with the size of a small custom microarray chip.

Table 1

Sorted white blood cells in S₀

Short name	Description	Number
B cells	CD19+ B-lymphocytes	6
Granulocytes	CD15+ granulocytes	8
Monocytes	CD14+ monocytes	5
NK	CD56+ Natural Killer (NK) cells	11
T cells (CD4+)1,2	CD3+CD4+ T-lymphocytes	8
T cells (CD8+)1,3	CD3+CD8+ T-lymphocytes	2
T cells (NKT)1	CD3+CD56+ natural killer	1
T cells (other)1	CD3+ T-lymphocytes	5

¹Considered as a member of the “pan-T-cell” group.

²Pan-T-cell further refined as also belonging to the “CD4+” group.

³Pan-T-cell further refined as also belonging to the “CD8+” group.

Figure 1

Clustering heatmap for external validation white blood cell data (S₀). Yellow = unmethylated (Y_hj = 0), black = partially methylated (Y_hj = 0.5), blue = methylated (Y_hj = 1).

Head and neck cancer

Our first target data set S₁ consisted of arrays applied to whole blood specimens collected in a random subset of individuals involved in an ongoing population-based case-control study[21] of head and neck cancer (HNSCC): 92 cases and 92 age and sex matched controls. The study was approved by Brown University IRB, protocol #0707992334. Blood was drawn at enrollment (prior to treatment in 85% of the cases). Mean age among the subjects arrayed in this study was 60 years, and there were 56 females and 128 males, consistent with the higher incidence of the disease in men. Thus, the covariate vector zconsisted of an indicator for case/control status, an indiator for male sex, and age (in decades) centered at the mean. The clustering heatmap in Figure2 depicts the raw DNA methylation data in S₁.

Figure 2

Clustering heatmap for target HNSCC data (S₁). Yellow = unmethylatedYij=0, black = partially methylatedYij=0.5, blue = methylatedYij=1. The annotation track above the heatmap indicates case-control status (orange = case, purple = control).

Simulations

We conducted extensive simulation studies in order to verify the finite-sample statistical properties of our proposed methodology. Simulation parameters were obtained from the HNSCC data set, and most simulations assumed no sources of biological bias (DNA methylation changes arising from processes not mediated by the profiled leukocytes, including shifts in distribution within cell types not profiled). In every simulation, we specified S₀ to consist of 5 B-cell samples, 10 granulocyte samples, 5 monocyte samples, 15 NK samples, 5 general “Pan-T” T-cell samples, 8 specific CD4+ T cell samples, and 2 specific CD8+ T cell samples. Estimates from the external validation set S₀, described above, were used for mean methylation profiles among WBC types, using the m = 100 most informative CpG sites.

Figure 3

Results of cell mixture reconstruction experiments validating prediction of individual profiles. Expected and observed percentages of each cell type are shown by color (red=100, white=0) and text. Median root-mean-square-error over 12 samples had a median value of 8.2%, ranging from 5.4% to 11.6%.

We specified n₁/2 cases and n₀/2 controls, n₀∈{100,200,500}. Among the controls, methylation profiles were generated by a white blood cell population of 7% B-cells, 62% granulocytes, 6% monocytes, 2% NK cells, and 13% were T-cells, of which 65% were CD4+ cells and 35% were CD8+ cells, and the remaining 5% were unspecified (and assumed to have mean methylation equal to that of the unsorted T-lymphocytes). Among cases, we specified one of the following scenarios: a 4% reduction in CD4+ cells, a 2% reduction in CD8+ cells, and an 8% increase in granulocytes (alternative with changes in both CD4+ and CD8+, “Strong Alternative I”); a 6% reduction in CD4+ cells, and an 8% increase in granulocytes (alternative with changes in CD4+ but not CD8+, “Strong Alternative II”); a weaker alternative with half the effects of Strong Alternative I (“Mixed Alternative” elaborated upon below); and two null scenarios with no changes in cell population, each with a different assumption about δ. Note that these changes reflect absolute changes in percentage points, not relative changes. Note also that these values were actually used to generate Dirichlet-distributed mixture weights for each simulated subject, with Dirichlet parameters equal to a precision parameter (100 corresponding to “precise” and 10 corresponding to“noisy”) times the mean weight described above. Residual effectsξi(0) for controls were set equal to 0.1 times estimated intercept estimateμ^1 obtained from the HNSCC data set, while residual effectsξi(1) for cases were set equal to 0.08 or 0.09 timesμ^1 plus multiples 10θ of the column ofU^corresponding to case. The constants of proportionality 0.1, 0.08, and 0.09 were chosen to correspond to assumed contributions of ξ to an overall methylation signature presumed to be dominated by profiled populations of white blood cells in specified proportions, with 0.08 used for the strong alternatives and 0.09 used for the Mixed Alternative. The constant 10 was used to amplify the scale of δ so that its effect could be detected in simulation; note thatU^was orthogonal to the white blood cell profiles, by construction. The multiplier θ = 0 was used for strong alternatives, and the “Strong Null” case (i.e. no methylation differences between cases and controls) while θ = 0.5 was used for the Mixed Alternative, and θ = 1 was used for the “Mixed Null” with case/control differences not mediated by cellular population differences. A simple normal error structure for e_0h and e_0i was specified, with no chip effects, but with variance equal to the sum of chip and residual variance estimated (individually for each CpG) for the HNSCC data. For each simulation, 50 bootstraps were used to estimate standard errors. 1000 simulations were run for each scenario.

Cell mixture experiment

As Figure3 suggests, accuracy is within 10%, and often less than 5%, with the largest errors occuring for granulocytes, as shown in Table2. Note that the sum of the individual observed predictions for each individual profile ranged from 98.9% to 102.7% (data not shown), even though the constraints of the projection do not explicitly constrain the sum to 100%; this provides additional evidence that the DNA methylation profile captures a great deal of information about cell mixtures.

Head and neck cancer

Table3 presents coefficient estimatesΓ^for case status, double-bootstrap bias estimates (estimates of bias arising from measurement error), as well as naive, single-bootstrap, and double-bootstrap standard error estimates. Each of these quantities is measured in percentage points (%). Estimates of bias arising from measurement error (i.e. substituting estimated quantities for known ones in a two-stage statistical procedure) were almost always less than half a percentage point, and for significant coefficient estimates, always towards the null. The proportion of CD4+ T-lymphocytes decreased in cases compared with controls, with a bias-corrected estimate of −10.4 percentage points and approximate 95% confidence interval (−13.1%,−3.3%); the proportion of NK cells decreased, with a bias-corrected estimate of -1.5 percentage points and 95% confidence interval (−2.2%,−0.75%); and the proportion of granulocytes increased, with a bias-corrected estimate of 7.6 percentage points and 95% confidence interval (4.2%,10.9%). There was also somewhat weaker evidence of an increase in CD8+ T-lymphocytes, with an estimate of 4.5 percentage points and 95% confidence interval (2.0%,7.0%). As reported in the complete set of results appearing in the Additional file1, the proportion of CD4+ T-lymphocytes decreased by 3.3 percentage points (−4.4%,−2.2%) per decade of age, while CD8+ T-lymphocytes increased by 2.0 percentage point (1.0%,3.0%) per decade. All other coefficients were insignificant.

For this analysis,R1,02 was estimated at 14.2%, whileR1,12 was estimated at 93.9%. Thus, a small but non-negligible proportion of total variation (systematic variation + unexplained biological heterogeneity + technical noise) appeared to be driven by changes in cell population between cases and controls and as a result of aging. Note that SS_e comprised 85% of total variation, so a substantial portion of variability in DNA methylation appeared to remain unexplained (presumably due, in large part, to technical noise). However, almost all of the systematic variation appeared to be explained by changes in cell population.

These results were consistent with previous studies, as HNSCC patients are known to display an absolute and relative increase in myeloid derived granulocytes[25] while also displaying an alteration in lymphoid T-cell homeostasis that leads to decreases in CD4+ T-cells[26,27]. In addition, the proportion of Treg cells (a subclass of CD4+ T cells) is known to decrease from infancy to adulthood[28].

The bias estimates obtained from the double-bootstrap procedure allow the correction of bias arising from measurement error. However, there is no statistical procedure for correcting the other possible sources of bias, those arising from changes in distribution among unprofiled cell types as well as non-immune-mediated methylation differences. The Additional file1 presents a detailed sensitivity analysis, from which we show that the magnitude of the resulting bias is likely to be small, less than a percentage point.

Ovarian cancer

Table4 presents results for case-control status, with the remaining results appearing in the Additional file1.R1,02 was estimated at 17.8%, whileR1,12 was estimated at 86.1%.

Compared with controls, cases showed significant increases in granulocytes and significant decreases in B cells, NK cells, and CD4+ T cells. Cases also showed marginally significant increases in monocytes. These results are consistent with previous literature, where it has been demonstrated that ovarian cancer patients experience decreases in B and T lymphocytes[29-31], increases in monocytes[29,30] and (somewhat equivocally) increases in eosinophil granulocytes[30]. Additionally, there were significant systematic decreases in CD4+ T cells with increasing age, with a gradient consistent in direction and somewhat consistent in magnitude with the corresponding effect found in the HNSCC data set. Though most of the CD8+ T cell coefficients for age were not significant, they were all positive, with gradient consistent in direction and somewhat consistent in magnitude with the corresponding effect found in the HNSCC data set. As reported in the Additional file1, no bisulfite conversion coefficient was significant, and all coefficients were of small magnitude (generally less than 1 percentage point per standard deviation).

Down syndrome

The only significant difference between cases and controls was in B cell distribution, with bias-corrected estimated decrease of 4.8%, 95% confidence interval (−6.2%, −3.5%). This result is consistent with known immune characteristics of Down Syndrome, including deficiencies in both B and T cells[32,33]. However, in the comparison between total leukocytes and T cells, all coefficients except B Cell and NK were highly significant, in directions consistent with comparison of a sample of purified T cells to a generic whole blood sample. In fact, an estimate of the cellular composition of the T cell samples can be obtained by a simple linear transformation of Γ estimates (adding intercept terms with the T cell coefficients); this operation produces values that are not significantly distinct from zero for all cell types except CD4+ and CD8+, whose bias-corrected estimates were, respectively, 75.9%, 95% confidence interval (67%,85%) and 8.6%, 95% confidence interval (0%,17%), consistent with the known distribution of these T cells. For the analysis of case vs. control within total leukocytes,R1,02 was estimated at 4.5%, whileR1,12 was estimated at 67.6%. For the analysis of total leukocyte vs. T cell with pooled cases and controls,R1,02 was estimated at 81.4%, whileR1,12 was estimated at 98.9%. The latter set of coefficients of determination indicate that a substantial portion of variation is explained by composition of leukocytes, which is the expected result for such an analysis.

Obesity in African Americans

Obese subjects had an estimated increase of 12 percentage points in granulocytes, bias-corrected 95% confidence interval (3.4%, 20%) and an estimated decrease of 4 percentage points in NK cells, bias-corrected 95% confidence interval (-7.7%,-0.9%). No significant differences were found for other blood cell types. Note that the specific immunological differences estimated by the method are consistent with known immunological perturbations associated with type II diabetes[9,10]. Complete results are provided in the Additional file1.

Additional analyses

We obtained the following unnormalized bias-corrected estimates: 69.0% CD4+, 95% CI (54%,84%), and 32.5% CD8+, 95%CI (19%,46%). This is consistent with known proportions of these specific cell types among T lymphocytes.

Results of simulations

Table5 presents results for n₁ = 200 with precise mixture weights (small within-status heterogeneity in distribution), while Table6 presents results for n₁ = 200 with noisy mixture weights (larger within-status heterogeneity). The tables show mean estimate, simulation standard deviation, median estimates for the three types of proposed standard errors, and proportion of p-values (obtained from z-scores constructed using the double-bootstrap standard error) falling below α = 0.05 and α = 0.01. In all cases, the bias in estimation was negligible. Both bootstrap procedures produced similar standard error estimates, which were close to the simulation standard deviation but often quite different from the naive standard error estimate. Under null scenarios, the rejection probabilities were tolerably close to their nominal values, and for alternatives, power could be quite high, even with this modest design. Results for the coefficients of determination are provided in the Additional file1. Scenarios with n₁∈{100,500} produced similar results, with simulation standard deviations and power adjusted accordingly, but still having practical utility.