Abstract
Background
There has been a longstanding need in biomedical research for a method that quantifies the normally mixed composition of leukocytes beyond what is possible by simple histological or flow cytometric assessments. The latter is restricted by the labile nature of protein epitopes, requirements for cell processing, and timely cell analysis. In a diverse array of diseases and following numerous immunetoxic exposures, leukocyte composition will critically inform the underlying immunobiology to most chronic medical conditions. Emerging research demonstrates that DNA methylation is responsible for cellular differentiation, and when measured in whole peripheral blood, serves to distinguish cancer cases from controls.
Results
Here we present a method, similar to regression calibration, for inferring changes in the distribution of white blood cells between different subpopulations (e.g. cases and controls) using DNA methylation signatures, in combination with a previously obtained external validation set consisting of signatures from purified leukocyte samples. We validate the fundamental idea in a cell mixture reconstruction experiment, then demonstrate our method on DNA methylation data sets from several studies, including data from a Head and Neck Squamous Cell Carcinoma (HNSCC) study and an ovarian cancer study. Our method produces results consistent with prior biological findings, thereby validating the approach.
Conclusions
Our method, in combination with an appropriate external validation set, promises new opportunities for largescale immunological studies of both disease states and noxious exposures.
Background
The biology of the development of any multisystem life form is fundamentally grounded in systematic cellular differentiation. This is essentially defined by lineage commitment of cells whose origin can be traced to a pluripotent progenitor and is marked by mitotically heritable epigenetic changes that reflect complex transcriptional programming of gene expression within the individual cell [13]. One such epigenetic mark is DNA methylation, which is tightly associated with alterations in the nucleosome DNA scaffold (and hence chromatin) that is responsible for coordination of gene expression in individual cells [13]. It is now appreciated that differentially methylated DNA regions (DMRs) distinguish cell lineages with high sensitivity and specificity [4] and considerable research is now underway to delineate precise DMRs that define and specify a particular cell lineage. The most developed understanding of epigenetic markers of lineage commitment to date is perhaps that of immune cell subclasses defined by populations of distinct circulating blood cells [5,6].
Pluripotent hematopoietic stem cells residing in the bone marrow continually give rise to the entire hierarchy of blood cell subclasses through a developmental process known as hematopoiesis. Leukocytes, commonly called white blood cells, are critical in the host response to pathogens and foreign antigens and are divided into two compartments, the myeloid lineage and lymphoid lineage (also called lymphocytes). The composition of leukocyte populations is well known to reflect disease states and toxicant exposures and can be altered by signaling cascades that prompt migration of whole classes of cells into or out of tissues. Several DMRs that serve as reliable biomarkers of individual human white blood cell types have already been identified [5,6]. Individual assays identifying cellspecific DMRs have proven useful for quantifying individual cell types in human tissues and peripheral blood. However, these assays are limited to detecting the relative proportion of one individual cell type compared with all others. On the other hand, simultaneous quantification of fluctuation in overall lymphocyte population composition can be accomplished only by using methods based on flow cytometry, which require large volumes of fresh blood and involve laborious antibody tagging. Hence, an approach that allows for the simultaneous quantification of the entire distribution of cell types, using an array of biomarkers based on generally available technology, would be considerably more informative, especially in studies of human disease and exposures.
In some instances, it is generally the overall balance of leukocyte subclasses in circulation or tissue that most prominently influences pathogenesis. For example, although incipient cancer cells are recognized and eliminated by cytotoxic Tcells (CTLs) and natural killer (NK) cells, tumorigenesis is also promoted by certain other inflammatory cells, including Blymphocytes, mast cells, neutrophils, regulatory Tcells (Tregs), and numerous others. All of these cells have been shown to promote angiogenesis, tumor cell proliferation, tissue invasion and metastasis [7,8]. Likewise, while higher levels of NK cells and CTLs circulating in the blood and residing in adipose tissues are associated with lower incidence of metabolic diseases such as type II diabetes [9], higher levels of M1 macrophages in adipose tissue can induce inflammation and insulin resistance [10]. These examples illustrate incredible potential for methods of quantifying the composition of lymphocyte populations to critically inform the underlying immunobiology of disease states as well as the immune response to almost all chronic medical conditions. In addition, they offer great potential for predicting therapeutic outcomes [11].
Here we employ the concept of DMRs as markers of immune cell identity using a high density methylation platform, and propose a set of analytical tools for estimating the proportions of immune cells in unfractionated whole blood that does not require fresh cells. The backbone of the approach is the DNA methylation signature of each of the principal immune components of whole blood (B cells, granulocytes, monocytes, NK cells, and T cells subsets). We essentially seek a form of regression calibration, where we consider a methylation signature to be a highdimensional multivariate surrogate for the distribution of white blood cells. In turn, this distribution is of interest for predicting or modeling disease states. As a surrogate, the DNA methylation signature is assumed to be a highly correlated, yet imperfect, measure of leukocyte distribution, and thus fits into the framework of measurement error models, where the use of a noisy surrogate marker to investigate an association with a disease outcome of interest results in biased estimates, unless internal or external validation data can be obtained to “calibrate” the model and correct the bias [12]. However, in this case, the problem is complicated by the extremely high dimension of the surrogate, so we propose an alternative to the traditional regressioncalibration procedure that circumvents these complications but still allows us to extract the desired biological information.
We note that since we began this work, a small number of authors have published similar deconvolution algorithms using gene expression data [1315]. The techniques are similar to the quadratic programming method we describe below in Methods for deconvolving a single sample, but none comprehensively addresses statistical properties or employs data from DNA methylation.
Methods
In this section we describe our proposed statistical methods, the data sets used to demonstrate their utility, and finally the design of simulation studies we have conducted to investigate statistical properties of our proposed algorithms.
Statistical methods
Let Y_{0h} be an m × 1 vector of methylation assay values, e.g. average beta values from an Infinium beadarray 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_{0} such types) indicated by a d_{0} × 1 covariate vector w_{h}. Here, h∈{1,…,n_{0}}, where n_{0} is the number of specimens and the m individual values correspond to CpG sites on a DNA methylation microarray, possibly preselected 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_{1}}, n_{1} is the number of target specimens, and z_{1i }is a d_{1}×1 covariate vector representing phenotypes or exposures corresponding to the subject, e.g. d_{1} = 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_{0} = {(Y_{01},w_{1}),…,(Y_{0n}_{0},w_{n}_{0})}, the set of data from “purified” cell samples effectively representing external validation or goldstandard data, and S_{1} = {(Y_{11},z_{1}),…,(Y_{1n}_{1},z_{n}_{1})}, representing surrogate data collected from a target population. To this end, we posit the following linear models:
where B_{0} and B_{1} are, respectively, m × d_{0} and m × d_{1} matrices and e_{0} and e_{1} are error vectors. For simplicity we assume a oneway ANOVA parameterization for w, though in the Additional file 1 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_{0} as μ_{1}, the m × 1 intercept. The error vectors e_{0} and e_{1} 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:
where Γ is a d_{0} × d_{1} 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_{0} = (b_{01},…,b_{0d}_{0}) explicitly in terms of its columns and writing
To impart a biological interpretation, we assume that the DNA assayed in S_{1} arises as a mixture of DNA from cell types profiled in S_{0}, with mixture coefficients whose population averages, conditional on z, are
Additional file 1. HousemanWBCBMCBioinformatics Supplement.pdf. Additional theoretical details, simulation descriptions and results, and additional figures and result tables [4352].
Format: PDF Size: 3.9MB Download file
This file can be viewed with: Adobe Acrobat Reader
where the m × 1 vector ξ^{(z)} represents cell types excluded from consideration among the purified samples in S_{0}, or else noncellspecific 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 Γ,
It is possible to assign interpretations to the components of variation in (3). Let
SS_{o} represents overall variability in Y_{1i}, i.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 estimates
Beyond bias due to measurement error, which is easily corrected using the doublebootstrap
procedure, there are additional sources of potential bias. For example, consider a
univariate z_{1i} representing case/control status, where
While the focus of this paper is analysis of population data, it is possible to use
S_{0} to predict distribution of leukocytes in a single sample having DNA methylation profile
Y^{∗}. Equating the intercept term of B_{1} in (1) with Y^{∗} and applying (2), we obtain mixing proportion estimates
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_{0} 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_{1}. 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_{1}. 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_{0} consisted of 46 white blood cell samples, deidentified 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 Table 1 and depicted by the clustering heatmap in Figure 1. Note that T lymphocytes that express CD4 or CD8 constitute over 95% of the T cell class, and that the panT cell type was further refined to CD4+, CD8+, and “other” PanT cells subtypes. In summary, the covariate vector w_{h}consisted of indicators for five cell types and another two indicators for CD4+ and CD8+ T cell subtypes. A generalization of the oneway ANOVA parameterization assumed above for w_{h}, described in the Additional file 1, was necessary to account for the ambiguous status of some PanT cells. For each CpG site, a linear mixed effects model with a random intercept for bead chip was used to estimate B_{0}; 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_{0}and B_{1}. From S_{0}, F statistics (described in the Additional file 1) 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 file 1, 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_{0}
Figure 1. Clustering heatmap for external validation white blood cell data (S_{0}). Yellow = unmethylated (Y_{hj }= 0), black = partially methylated (Y_{hj }= 0.5), blue = methylated (Y_{hj }= 1).
Cell mixture experiment
Proof of the utility of the proposed methods in predicting leukocyte distributions for individual samples requires extensive, detailed reconstruction experiments beyond the scope of the present paper. However, to provide evidence that such experiments are worthwhile and show promise of positive results, we conducted a simple experiment involving six known mixtures of monocytes and B cells and six known mixtures of granulocytes and T cells. The results of this experiment are described below in Results.
Head and neck cancer
Our first target data set S_{1} consisted of arrays applied to whole blood specimens collected in a random subset of individuals involved in an ongoing populationbased casecontrol 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 Figure 2 depicts the raw DNA methylation data in S_{1}.
Figure 2. Clustering heatmap for target HNSCC data (S_{1}). Yellow = unmethylated
Ovarian cancer
We next applied our method to an ovarian cancer data set [22]. DNA methylation data for blood samples are available from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/ webcite, Accession number GSE19711). We used only those cases having blood drawn pretreatment. After removing 4 arrays with a preponderance of missing values, the data set consisted of 272 controls and 129 cases having blood drawn prior to treatment. A clustering heatmap displaying the DNA methylation data appears in the Additional file 1. In this analysis, zconsisted of casecontrol status, age (categorized in 5year increments), and 2 bisulfite conversion efficiency measures.
Down syndrome
We also applied our method to a trisomy 21 (Down syndrome) data set [23] consisting of 29 total peripheral blood leukocyte samples from Down syndrome cases and 21 controls, as well as 6 T cell samples from cases and 4 T cell samples from controls (GEO Accession number GSE25395). Because of the potential for bias induced by copy number amplification, we excluded 4 CpG sites on Chromosome 21, resulting in m = 96 CpG sites used for analysis. A clustering heatmap displaying the DNA methylation data appears in the Additional file 1. In one analysis, we compared cases and controls using the total leukocyte samples only, and in another we compared total leukocytes to T cells, pooling cases and controls. The Additional file 1. presents coefficient estimates.
Obesity in African Americans
Finally, we applied our method to an obesity data set [24] consisting of 7 lean AfricanAmericans and 7 Obese AfricanAmericans (GEO Accession number GSE25301). A clustering heatmap displaying the DNA methylation data appears in the Additional file 1. In this analysis, zconsisted of obesity status.
Additional analyses
If the subject population for which z = 0 is sufficiently homogeneous with respect to blood cell distribution to admit sensible
characterization of that distribution, then it is possible to recover estimates from
Simulations
We conducted extensive simulation studies in order to verify the finitesample 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_{0} to consist of 5 Bcell samples, 10 granulocyte samples, 5 monocyte samples, 15 NK samples, 5 general “PanT” Tcell samples, 8 specific CD4+ T cell samples, and 2 specific CD8+ T cell samples. Estimates from the external validation set S_{0}, 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 rootmeansquareerror over 12 samples had a median value of 8.2%, ranging from 5.4% to 11.6%.
We specified n_{1}/2 cases and n_{0}/2 controls, n_{0}∈{100,200,500}. Among the controls, methylation profiles were generated by a white
blood cell population of 7% Bcells, 62% granulocytes, 6% monocytes, 2% NK cells,
and 13% were Tcells, 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 Tlymphocytes). 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 Dirichletdistributed
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
Results
In this section we report the results of the data analyses described above in Implementation, as well as the results of our simulation experiments.
Cell mixture experiment
As Figure 3 suggests, accuracy is within 10%, and often less than 5%, with the largest errors occuring for granulocytes, as shown in Table 2. 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.
Table 2. Summary statistics for errors in cell mixture reconstruction results^{*}
Head and neck cancer
Table
3 presents coefficient estimates
For this analysis,
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 Tcell homeostasis that leads to decreases in CD4+ Tcells [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 doublebootstrap 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 nonimmunemediated methylation differences. The Additional file 1 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
Table
4 presents results for casecontrol status, with the remaining results appearing in
the Additional file
1.
Table 3. Estimates for HNSCC analysis (case vs. control)
Table 4. Estimates for ovarian cancer analysis (case vs. control)
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 [2931], 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 file 1, 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 biascorrected 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 biascorrected 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,
Obesity in African Americans
Obese subjects had an estimated increase of 12 percentage points in granulocytes, biascorrected 95% confidence interval (3.4%, 20%) and an estimated decrease of 4 percentage points in NK cells, biascorrected 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 file 1.
Additional analyses
We obtained the following unnormalized biascorrected 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
Table 5 presents results for n_{1} = 200 with precise mixture weights (small withinstatus heterogeneity in distribution), while Table 6 presents results for n_{1} = 200 with noisy mixture weights (larger withinstatus heterogeneity). The tables show mean estimate, simulation standard deviation, median estimates for the three types of proposed standard errors, and proportion of pvalues (obtained from zscores constructed using the doublebootstrap 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 file 1. Scenarios with n_{1}∈{100,500} produced similar results, with simulation standard deviations and power adjusted accordingly, but still having practical utility.
Discussion
In this paper, we employ the concept of DMRs as markers of immune cell identity using a high density methylation platform, and propose a set of analytical tools for estimating the proportions of immune cells in unfractionated whole blood. The backbone of the approach is the DNA methylation signature of each of the principal immune components of whole blood (B cells, granulocytes, monocytes, NK cells, and T cells subsets). The examples we have provided above serve to illustrate that our proposed methodology produces parameter estimates consistent with the literature, thus validating its utility.
Our proposed method resembles regression calibration, where we consider a methylation signature to be a highdimensional multivariate surrogate for the distribution of white blood cells. In turn, this distribution is of interest for predicting or modeling disease states. As a surrogate, the DNA methylation signature is assumed to be a highly correlated, yet imperfect, measure of leukocyte distribution, and thus fits into the framework of measurement error models, where the use of a noisy surrogate marker to investigate an association with a disease outcome of interest results in biased estimates, unless internal or external validation data can be obtained to “calibrate” the model and correct the bias [12]. However, in this case, the problem is complicated by the extremely high dimension of the surrogate. Measurement error problems are typically formulated as a set of relationships between z, the disease outcome (e.g. case/control status), ω, the gold standard (e.g. leukocyte distribution), and Y, the surrogate (e.g. DNA methylation). Of interest is E(z  ω), which may be difficult to estimate due to the cost or logistical complications involved in obtaining ω in a large number of samples. Typically, it is possible to collect sufficient data for modeling E(z  Y), which provides information about E(z  ω) through the (often imperfect) association E(Y  ω), which is inferred from an external validation sample [12,34]. Unfortunately, the highdimensional nature of Y renders E(z  Y) difficult to formulate. While multivariate methods of measurement error correction exist, even in a highdimensional context [35], they require an explicit specification of E(zY), requiring a large number of parameters even for a main effects regression model, and many more in order to account for interactions. This becomes unwieldy when each component of Y contributes a small amount of information about z, and both dimensionreduction strategies and constrained regression strategies entail substantial loss of information and may be extremely computationally intensive. Existing measurement error formulations [34,35] would have required us to specify a logistic regression model for case/control status, conditional on DNA methylation signature, a computationally difficult task that would have extreme vulnerability to model misspecification. On the other hand, our method requires specification of E(Y  z), which is natural and straightforward. Note that in some treatments of regression calibration, E(ω  Y) is used as a surrogate for ω in regression models for z[12]; our treatment essentially assumes a linear form for E(Y  ω) and effectively obtains E(ω  Y) by projecting Y onto the column space of resulting matrix. We note that it is possible using existing methods to qualitatively describe immune response contributions to DNA methylation. This is typically done by conducting a pathway analysis along the lines of one of the methods described in [36], the best option of which is Gene Set Enrichment Analysis (GSEA) [37]. For example, Teschendorff et al. (2009) [22] use GSEA to qualitatively motivate an immunological explanation. However, these methods do not directly quantify the immunological contribution.
An important consideration in the measurement error literature is that of transportability of model parameters [38]. In our setting, an important consideration is whether the methylation profiles obtained from the purified blood cells used to assemble S_{0} would be representative of the white blood cells measured within S_{1}. Because of the biological assumptions inherent in the DMR literature and underlying current understanding of hematopoeisis and lineage commitment, this assumption is reasonable, provided our method is used to detect abnormal mixtures of normal white blood cells. However, methylation abnormalities in the white blood cells themselves constitute a form of noncell mediated alteration (in the sense of the term we have been using), and contribute to bias in our methods, as described briefly above and in detail in the Additional file 1.
Note that our formulation respects the study design (DNA methylation assay data collected after sampling from phenotype groups). An alternative strategy outside the measurement error literature but within the larger missingdata literature might have been the use of an ExpectationMaximization (EM) algorithm to integrate over the missing data ω[39]. However, by design, the distribution of ω varied substantially between the data sets S_{0} and S_{1}, severely complicating the approach; notably, an would be the introduction of feedback from S_{1} to S_{0}, contaminating the goldstandard status of S_{0}. An alternative, might be the use of an empirical Bayes procedure, reminiscent of existing mixturemodel approaches [40]. However, difficulty in specifying the distribution of “remainder terms” (denoted as ξabove) render this approach untenable, and in simulations (not presented), attempts to impute ωamong S_{1} samples using parameters obtained from S_{0} samples resulted in extremely biased estimates of ω.
The most significant aspect of the current study is our development of a method for inferring changes in the distribution of white blood cell types between different human populations (e.g. cases and controls) using DNA methylation signatures; an approach guided by an external validation set consisting of methylation profiles from purified white blood cell components. DNA methylation in peripheral blood is a potentially powerful new biomarker for clinical and epidemiological investigation. By example, numerous studies have now attempted to distinguish cancer cases from controls using whole peripheral blood assayed via DNA methylation arrays, including ovarian [22], bladder [41], and pancreatic [42] cancers. While these studies have demonstrated good to excellent discrimination of cases from controls, sound evidence for a biological mechanism has been elusive. Presumably, disease associated alterations in blood methylation have several etiological components driven by inherent genetic, environmental and disease specific factors. Given the known developmental associated differences in DNA methylation among specific blood cell types, changes in the distributions of blood cell types alone could account for disease associated DNA methylation. While numerous authors provide a qualitative discussion that includes the possibility of immunerelated DNA methylation differences (e.g. [22]), none to date has specifically quantified the contribution from immune response. On the other hand, the many diverse types of immune cells in blood make this issue highly complex and problematic to tackle using single cell type assays. Therefore, it is crucial to the development of this new avenue of biomarker research to delineate effects due to the immune cell distribution itself from other “non cell type” alterations in DNA methylation. We term the differences among human populations attributed to cell distributions to be “immunologically mediated”. Our solution to partition this component of variation in methylation from other determinants are multivariate analytic tools including regression coefficients and associated inference, as well as coefficients of determination measures. Taken together these provide a means for evaluating whether the observed DNA methylation differences are due to an immunologically mediated response.
In our Additional file 1 we provide a detailed analysis of potential sources of bias in our analysis. One obvious biological source of bias is age of the subjects contributing cells for validation. At certain CpG loci, DNA methylation is known to change with age [43], especially in T cells [44]. In the Additional file 1 we demonstrate that any agerelated associations with DNA methylation in our top 100 CpGs were too weak to be detected with the current validation sample, and thus unlikely to bias the results of our analyses (notably age coefficients provided for the HNSCC example). However, we remark that with larger sample sizes, adjustments for age can be incorporated with an appropriate additional term in the linear model (1) for Y_{0h}.
Similar methods based on mRNA have been employed [1315]. The statistical principles described in this article would apply, wholesale, to mRNA expression profiles, but with two cautionary statements. The first is mathematical: mRNA is typically analyzed on a logarithmic scale, yet the assumptions of the proposed methodology involve linearity on an arithmetic scale, since the mixing coefficients are assumed to act linearly on absolute numbers of nucleic acid molecules; thus, the proposed methods would require analysis of untransformed fluorescence intensities, whose skewed distributions would result in numerical instabilities. The second is biological: there is no necessarily linear relationship between cell number and mRNA copies, since proteins may be translated as a consequence of an initial burst of mRNA transcription upon cellular development, after which significant mRNA degradation is possible. In contrast, one would expect the average beta value provided by Illumina beadarray products (and similar quantities) to scale in proportion to the actual fraction of methylated nucleic acids; in addition, an assumption of two DNA molecules per cell seems biologically reasonable. In the Additional file 1 we provide an example of an application of our methods using mRNA data.
Going forward there are two issues that require further experimental and analytical refinement. First, although the current studies suggest group level comparisons of blood cell DNA methylation can reveal important immune alterations, it will be important to provide methods for individual level immune cell profiling, since clinical and detailed analytical epidemiologic applications that examine individual risk factor information will be the subject of future studies. As we have demonstrated above, individual immune profiles are theoretically achievable but will require extensive validation, with a wide array of mixture combinations, before gaining widespread acceptance. Secondly, there is intense interest in minor immune cell fractions and their role in disease, though the signal strength of cell types comprising < 5% of the total white cell compartment may be difficult to quantitate. Examples of such cell types include the regulatory Tcell or NK cell fractions, which are implicated in autoimmune and malignant diseases. Optimization of platforms for technical sensitivity to minor subtypes combined with statistical optimization of signature recognition are needed to enhance the approach for testing highly targeted immune hypotheses.
Conclusions
The method we present here has potentially far reaching implications for rapid, simple and complete assessment of the composition of human white blood cell populations, i.e. the immune profile. Currently, assessment of the cellular composition of peripheral blood cannot be accomplished without the use of freshly drawn venous blood that is immediately prepared in a specially equipped laboratory. A complete assessment of the entire immune profile requires extensive flow cytometric measurements based on protein epitopes on leukocyte membranes that distinguishes subtypes of immune cells that are either too rare or too similar in appearance to be distinguished using simple microscopic approaches. In particular, flow cytometry is limited by the following: (i) cells must be separated, requiring large volumes of fresh cells; (ii) detection can be accomplished only by the fluorescent antibody tags available, which require expensive technology to read; (iii) the outer cell membrane must be intact, mandating limited utility in many instances (particularly in research). In contrast, our method requires the application of these laborintensive or expensive steps only in the construction of the validation set S_{0}, which need only be developed once. Once S_{0} is available, subsequent interrogation is based on the chemically stable CpG methylation of DNA; thus our method obviates the need for fresh blood and the preservation of labile protein epitopes. It is also able to simultaneously assess all of the individual components of the peripheral blood using a highly multiplexed molecular platform and is thus very straightforward logistically. Furthermore, the statistical methodology presented here can be implemented easily with the instrumental output of the methylation arrays, which simplifies the interpretation of the immune profile data from the operators point of view. This method can be immediately deployed in a research framework to cost effectively assess human immune profiles (in fresh or archival samples), exploring their potential as biomarkers, and addressing key questions regarding disease pathogenesis. Furthermore, our approach is readily suited for rapid translation to a broad base of clinical applications such as disease monitoring, diagnosis, prognosis, and response to therapy.
Our approach makes research on biobanked specimens possible, now making a vast array of prospective studies that could not otherwise be done, possible. Software and sample data are provided in Additional file 2.
Additional file 2. HousemanWBCBMCBioinformaticsSoftwarev2. Sample R software (compressed).
Format: ZIP Size: 969KB Download file
Abbreviations
CTL: Cytotoxic Tcells; CpG: Cytosinephosphateguanine; DMR, Differentially methylated region; HNSCC: Head and neck squamous cell carcinoma; NK: Matural killer.
Competing interests
A patent is pending on the work contained in this article. The authors have no other competing interests.
Authors’ contributions
EAH conceived of the statistical model, developed the algorithms, conducted the simulations, applied the methods to proprietary and publicly available data sets and authored major parts of the manuscript. WPA conducted the laboratory experiments and authored parts of the manuscript. JKW and KTK conceived of the laboratory experiments and provided grant support for the research. DCK provided indispensible feedback on statistical methodology. BCC, CJM, and HHN provided indispensible feedback on scientific issues and interpretation. All authors read and approved the final manuscript.
Acknowledgements
This work was funded by NIH grants CA100679, CA126939, CA121147, and CA078609.
References

Natoli G: Maintaining cell identity through global control of genomic organization.

Ji H, Ehrlich LI, Seita J, Murakami P, Doi A, Lindau P, Lee H, Aryee MJ, Irizarry RA, Kim K, Rossi DJ, Inlay MA, Serwold T, Karsunky H, Ho L, Daley GQ, Weissman IL, Feinberg AP: Comprehensive methylome map of lineage commitment from haematopoietic progenitors.

Khavari DA, Sen GL, Rinn JL: DNA methylation and epigenetic control of cellular differentiation.

Baron U, Turbachova I, Hellwag A, Eckhardt F, Berlin K, Hoffmuller U, Gardina P, Olek S: DNA methylation analysis as a tool for cell typing.
Epigenetics 2006, 1:5560. PubMed Abstract  Publisher Full Text

Wieczorek G, Asemissen A, Model F, Turbachova I, Floess S, Liebenberg V, Baron U, Stauch D, Kotsch K, Pratschke J, Hamann A, Loddenkemper C, Stein H, Volk HD, Hoffmuller U, Grutzkau A, Mustea A, Huehn J, Scheibenbogen C, Olek S: Quantitative DNA methylation analysis of FOXP3 as a new method for counting regulatory T cells in peripheral blood and solid tissue.
Cancer Res 2009, 69(2):599608. PubMed Abstract  Publisher Full Text

Sehouli J, Loddenkemper C, Cornu T, Schwachula T, Hoffmuller U, Grutzkau A, Lohneis P, Dickhaus T, Grone J, Kruschewski M, Mustea A, Turbachova I, Baron U, Olek S: Epigenetic quantification of tumorinfiltrating Tlymphocytes.
Epigenetics 2011, 6(2):236246. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hanahan D, Weinberg RA: Hallmarks of cancer: the next generation.
Cell 2011, 144(5):64674. PubMed Abstract  Publisher Full Text

OstrandRosenberg S: Immune surveillance: a balance between protumor and antitumor immunity.
Curr Opin Genet Dev 2008, 18:1118. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lynch LA, O’Connell JM, Kwasnik AK, Cawood TJ, O’Farrelly C, O’Shea DB: Are natural killer cells protecting the metabolically healthy obese patient?
Obesity 2009, 17(3):601605. PubMed Abstract  Publisher Full Text

Anderson EK, Gutierrez DA, Hasty AH: Adipose tissue recruitment of leukocytes.

Chua W, Charles KA, Baracos VE, Clarke SJ: Neutrophil/lymphocyte ratio predicts chemotherapy outcomes in patients with advanced colorectal cancer.
Brit J Cancer 2011, 104:12881295. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Carroll RJ, Ruppert D, Stefanski LA: Measurement Error in Nonlinear Models. Boca Raton, Florida: Chapman & Hall; 2006.

Gaujoux R, Seoighe C: Semisupervised Nonnegative Matrix Factorization for gene expression deconvolution: a case study.
Infect Genet Evol 2011.
[10.1016/j.meegid.2011.08.014]

Gong T, Hartmann N, Kohane IS, Brinkmann V, Staedtler F, Letzkus M, Bongiovanni S, Szustakowski JD: Optimal deconvolution of transcriptional profiling data using quadratic programming with application to complex clinical blood samples.
PLoS One 2011, 6:e27156. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

ShenOrr SS, Tibshirani R, Khatri P, Bodian DL, Staedtler F, Perry NM, Hastie T, Sarwal MM, Davis MM, Butte AJ: Cell typespecific gene expression differences in complex tissues.

Wang SC, Petronis A: DNA Methylation Microarrays: Experimental Design and Statistical Analysis. Boca Raton, Florida: Chapman & Hall; 2008.

Smyth GK: Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.

Leek JT, Storey JD: Capturing heterogeneity in gene expression studies by surrogate variable analysis.
PLoS Genet 2007, 3:17241735. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Teschendorff AE, Zhuang J, Widschwendte rM: Independent surrogate variable analysis to deconvolve confounding factors in largescale microarray profiling studies.
Bioinformatics 2011, 27(11):14961505. PubMed Abstract  Publisher Full Text

Goldfarb D, Idnani A: A numerically stable dual method for solving strictly convex quadratic programs.
Math Prog 1983, 27:133. Publisher Full Text

Peters ES, McClean MD, Liu M, Eisen EA, Mueller N, Kelsey KT: The ADH1C polymorphism modifies the risk of squamous cell carcinoma of the head and neck associated with alcohol and tobacco use.
Cancer Epidemiol Biomarkers Prev 2005, 14(2):476482. PubMed Abstract  Publisher Full Text

Teschendorff AE, Menon U, GentryMaharaj A, Ramus SJ, Gayther SA, Apostolidou S, Jones A, Lechner M, Beck S, Jacobs IJ, Widschwendter M: An epigenetic signature in peripheral blood predicts active ovarian cancer.
PLoS ONE 2009, 4(12):e8274. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kerkel K, Schupf N, Hatta K, Pang D, Salas M, Kratz A, Minden M, Murty V, Zigman WB, Mayeux RP, Jenkins EC, Torkamani A, Schork NJ, Silverman W, Croy BA, Tycko B: Altered DNA methylation in leukocytes with trisomy 21.
PLoS Genet 2010, 6(11):e1001212. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wang X, Zhu H, Snieder H, Su S, Munn D, Harshfield G, Maria BL, Dong Y, Treiber F, Gutin B, Shi H: Obesity related methylation changes in DNA of peripheral blood leukocytes.
BMC Med 2010, 8:87. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Trellakis S, Bruderek K, Dumitru CA, Gholaman H, Gu X, Bankfalvi A, Scherag A, Hutte J, Dominas N, Lehnerdt GF, Hoffmann TK, Lang S, Brandau S: Polymorphonuclear granulocytes in human head and neck cancer: Enhanced inflammatory activity, modulation by cancer cells and expansion in advanced disease.
Int J Cancer 2011.
[10.1002/ijc.25892]

Kuss I, Hathaway B, Ferris RL, Gooding W, Whiteside TL: Decreased absolute counts of T lymphocyte subsets and their relation to disease in squamous cell carcinoma of the head and neck.
Clin Cancer Res 2004, 10(11):37553762. PubMed Abstract  Publisher Full Text

Kuss I, Hathaway B, Ferris RL, Gooding W, Whiteside TL: Imbalance in absolute counts of T lymphocyte subsets in patients with head and neck cancer and its relation to disease.
Adv Otorhinolaryngol 2005, 62:161172. PubMed Abstract  Publisher Full Text

Mold JE, Venkatasubrahmanyam S, Burt TD, Michaelsson J, Rivera JM, Galkina SA, Weinberg K, Stoddart CA, McCune JM: Fetal and adult hematopoietic stem cells give rise to distinct T cell lineages in humans.
Science 2010, 330(6011):16951699. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

den Ouden M, Ubachs JMH, Stoot JEGM, van Wersch JWJ: Whole blood cell counts and leucocyte differentials in patients with benign or malignant ovarian tumours.
Eur J Obstet Gynecol Reprod Biol 1997, 72:7377. PubMed Abstract  Publisher Full Text

Bishara S, Griffin M, Cargill A, Bali A, Gore ME, Kaye SB, Shepherd JH, Van Trappen PO: Pretreatment white blood cell subtypes as prognostic indicators in ovarian cancer.

Cho H, Hur HW, Kim SW, Kim SH, Kim JH, Kim YT, Lee K: Pretreatment neutrophil to lymphocyte ratio is elevated in epithelial ovarian cancer and predicts survival after treatment.
Cancer Immunol Immunother 2009, 58:1523. PubMed Abstract  Publisher Full Text

Verstegen RH, Kusters MA, Gemen EF, De Vries E: Down syndrome Blymphocyte subpopulations, intrinsic defect or decreased Tlymphocyte help.
Pediatr Res 2010, 67:563569. PubMed Abstract  Publisher Full Text

Ram G, Chinen J: Infections and immunodeficiency in Down syndrome.
Clin Exp Immunol 2011, 164:916. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Thurston SW, Spiegelman D, Ruppert D: Equivalence of regression calibration methods for main study/external validation study designs.
J Stat Plan Inf 2003, 113:527534. Publisher Full Text

Li B, Yin X: On surrogate dimension reduction for measurement error regression: an invariance law.
Ann Stat 2007, 35(5):21432172. Publisher Full Text

Goeman J, Buhlmann P: Analyzing gene expression data in terms of gene sets: methodological issues.
Bioinformatics 2007, 23:980987. PubMed Abstract  Publisher Full Text

Subramanian A, Tamayo P, Mootha V, Mukherjeed S, Ebert B, Gillette M, Paulovich A, Pomeroy S, Golub T, Lander E, Mesirov J: Gene set enrichment analysis: a knowledgebased approach for interpreting genomewide expression profiles.
Proc Natl Acad Sci USA 2005, 102:1554515550. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Carroll RJ, Galindo CD: Measurement error, biases, and the validation of complex models for blood lead levels in children.
Env Health Persp 1998, 106:15351539. Publisher Full Text

Little RJA, Rubin DB: Statistical Analysis with Missing Data. Hoboken, NJ: Wiley; 2002.

Koestler DC, Marsit CJ, Christensen BC, Karagas MR, Bueno R, Sugarbaker DJ, Kelsey KT, Houseman EA: Semisupervised recursively partitioned mixture models for identifying cancer subtypes.
Bioinformatics 2010, 26(30):25782585. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Marsit CJ, Koestler DC, Christensen BC, Karagas MR, Houseman EA, Kelsey KT: DNA methylation array analysis identifies profiles of bloodderived DNA methylation associated with bladder cancer.
J Clin Oncol 2011, 29(9):11331139. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pedersen KS, Bamlet WR, Oberg AL, de Andrade M, Matsumoto ME, Tang H, Thibodeau SN, Petersen GM, Wang L: Leukocyte DNA methylation signature differentiates pancreatic cancer patients from healthy controls.
PLoS ONE 2011, 6(3):e18223. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bocklandt S, Lin W, Sehl ME, Sanchez FJ, Sinsheimer JS, Horvath S, Vilain E: Epigenetic predictor of age.
PLoS One 2011, 6(6):e14821. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chu M, Siegmund KD, Hao QL, Crooks GM, Tavare S, Shibata D: Inferring relative numbers of human leucocyte genome replications.

Doi A, Park IH, Wen B, Murakami P, Aryee MJ: Differential methylation of tissue and cancerspecific CpG island shores distinguishes human induced pluripotent stem cells, embryonic stem cells and fibroblasts.
Nat Genet 2009, 41:13501353. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Houseman EA, Christensen BC, Yeh RF, Marsit CJ, Karagas MR: Modelbased clustering of dna methylation array data: a recursivepartitioning algorithm for highdimensional data arising as a mixture of beta distributions.
BMC Bioinf 2008, 9:365. BioMed Central Full Text

Alberts B, Johnson A, Lewis J, Raff M, Roberts K: Molecular Biology of the Cell. New York, NY: Taylor & Francis; 2008.

Showe MK, Vachani A, Kossenkov AV, Yousef M, Nichols C: Gene expression profiles in peripheral blood mononuclear cells can distinguish patients with nonsmall cell lung cancer from patients with nonmalignant lung disease.
Cancer Res 2009, 69:92029210. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kossenkov AV, Vachani A, Chang C, Nichols C, Billouin S: Resection of nonsmall cell lung cancers reverses tumorinduced gene expression changes in the peripheral immune system.
Clin Cancer Res 2011, 17:58675877. PubMed Abstract  Publisher Full Text

Watkins NA, Gusnanto A, de Bono B, De S, MirandaSaavedra D: A haematlas: characterizing gene expression in differentiated human blood cells.
Blood 2009, 113:e1e9. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ginns LC, Goldenheim PD, Miller LG, Burton RC, Gillick L: Tlymphocyte subsets in smoking and lung cancer: analysis of monoclonal antibodies and flow cytometry.

Mazzoccoli G, Balzanelli M, Giuliani A, De Cata A, La Viola M: Lymphocyte subpopulations anomalies in lung cancer patients and relationship to the stage of disease.
In Vivo 1999, 13:205209. PubMed Abstract