Measuring cell-type specific differential methylation in human brain tissue

Estimation of mixture proportions

We measured DNA methylation profiles for dorsolateral prefrontal cortex (DLPFC), hippocampal formation (HF), and superior temporal gyrus (STG) samples dissected from frozen brains of normal individuals using the comprehensive high-throughput arrays for relative methylation (CHARM) technique [23]. We also labeled and separated neuronal nuclei in a subset of samples using a neuron-specific antibody (NeuN) and fluorescence-activated cell sorting (FACS) [24, 25]. Neuronal (NeuN+) and non-neuronal (NeuN-) fractions from DLPFC, HF, and STG were collected for downstream processing and methylation analysis with CHARM (Additional File 1, Figure S1).

To illustrate the downstream effects of the cell population confounding problem, and focusing on two brain regions for clarity, we examined a genomic region for which: (1) no difference was observed between DLPFC and HF in either neuronal or glial fractions; and (2) a difference was observed between neuronal and glial nuclei within each brain region (Figure 1a). Note that a strong methylation difference between brain regions is observed between the non-cell-sorted brain samples. This must be a false-positive and, as we demonstrate below, must be due to differences in cell-type composition between the brain regions.

We modified a statistical method originally developed to estimate cell populations in blood [19] to calculate neuronal and glial proportions for each of our unsorted samples, adapting it to use a constrained linear optimization model (Figure 1b, see overview in Additional File 1, Figure S2a). We confirmed that our approach effectively estimated these cell proportions using a mixture experiment with an independent set of samples (Additional File 1, Figure S2b). To demonstrate that the false-positive results of Figure 1 are due to difference in cell-type distribution, we mathematically reconstructed the unsorted sample methylation profile using the pure neuronal and glial profiles and their estimated frequencies and predicted this result (Additional File 1, Figure S2c).

While the above results rely on having neuronal and glial methylation signals for each brain region, we performed additional analyses to determine whether accurate estimates of neuronal and glial proportions in unsorted samples from a brain region could be obtained using selected data from another brain region. Figure 1c shows the accuracy of estimates obtained from such 'universal' data, compared to estimates based on sorted data from each individual brain region. We also accurately reproduce the cell proportion estimates from our mixture experiment (Additional File 1, Figure S2d, see Materials and Methods for additional details of how this analysis was performed). Our results indicate that accurate estimates could be obtained for a new brain region without the need to sort samples from that region.

Generative model of methylation signal

in which the parameters ${\beta }_2$ and ${\beta }_3$ represent the quantities we are interested in measuring, that is, the differences in neurons and glia, respectively, between brain regions. We refer to this model as M2. Fitting this linear model by least squares and obtaining estimates for millions of genomic locations is computationally feasible. (Fitting the model for 4 million probes took about 5 seconds on our laptop).

This statistical framework also exposes the problem with existing naïve approaches to assess DNA methylation signatures in mixed samples. To date, most published analyses ignore cell composition [26–30] and look for associations in a way equivalent to fitting a simple linear regression model $Y_i={\alpha }_0+{\alpha }_1{X}_i+{\epsilon }_i$ (where the t-test is derived from the $X_i=0\mathsf{\text{or}}1$). We refer to this model as M1. In M1, the parameter ${\alpha }_1$ represents a combination of the methylation differences in neurons and glia in which it is impossible to deconvolve cell-type-specific contributions. Furthermore, we can mathematically demonstrate that the least squares estimate of ${\alpha }_1$ will be biased by differences in cell-type frequency under the null hypothesis of no difference in methylation between brain regions (Figure 2a, see Methods Section). Similarly, a naïve model suggested by Guintivano et al. [20] that incorporates cell-type proportions $Y_i={\gamma }_0+{\gamma }_1{X}_i+{\gamma }_2{\pi }_i+{\epsilon }_i$ (we refer to this as model M3) will lead to biased results as well, and to decreased power to detect methylation differences (Additional File 1, Figure S3). We also note that even the superior methods show a small amount of bias (boxplot not centered at 0), which can be explained by slightly inaccurate mixture estimates (see Materials and Methods).

To test the utility of our model, we confirmed our theoretical results with experimental data. First, we obtained estimates of significant neuron-specific methylation differences between DLPFC and HF using sorted brain samples (gold standard, FDR <0.05, Additional File 1, Table S1). We then used the unsorted brain data to calculate the parameters representing the differences in brain-region methylation using models M1 (total methylation difference, ${\alpha }_1$) and M2(neuron-specific methylation difference, ${\beta }_2$). Figure 2b shows that we can estimate neuron-specific methylation differences more accurately with model M2. Therefore, we can assess neuron-specific methylation differences between DLPFC and HF using whole tissue after estimating cell proportions.

Using the sorted samples, we did not find statistically significant DMRs in the non-neuronal fraction, which highlights the importance of isolating a neuronal signal from total methylation values. The result is in agreement with recently published literature suggesting that glia cells, contained in the NeuN- fraction, have less diverse transcription patterns across brain regions than neurons [31], the latter having a distinct DNA-methylation signature [10]. Interestingly, proteins involved in modifying chromatin were found among the brain-region neuronal DMRs, supporting the role of epigenetic mechanisms in neuronal function and synaptic plasticity [32]. For example, neuron-specific methylation of the histone methyltransferase SETD3, which methylates histone H3 at lysine 36, was lower in HF than in DLPFC, and histone deacetylase HDAC4 shows hypomethylation in DLPFC. Other genes involved in neural differentiation include JAG1, TTL1, NPAS4, CUX-2, DOCK2, NGEF, OLFM1, SATB2, and GIT2.

Application to Illumina Infinium HumanMethyation450 Dataset

While the CHARM platform has many advantages for studying methylation patterns due to the high density and location of probes, the assay requires restriction-enzyme digestion and lacks single-base resolution. The Illumina Infinium HumanMethylation450 (450K) array has emerged as an affordable alternative to obtain reliable quantitative measurements of methylation. To demonstrate the performance of our method on data from the 450K array, we used data accessible at NCBI GEO database (Guintivano et al. [20], accession GSE41826), consisting of 77 normal samples from prefrontal cortex, of which 29 were sorted into neuronal and glial fractions, nine were mixtures of neurons and glia of known proportions, and 10 were unsorted, whole-tissue samples. We first applied our method to obtain accurate cell-fraction estimates on the known mixture samples (Additional File 1, Figure S4a). Using these cell-fraction estimates and the pure neuronal and glial profiles, we mathematically reconstructed the methylation profile for the mixture samples in a set of genomic regions and compared these results to the actual observed methylation for these samples (Additional File 1, Figure S4b). The cell proportion calculations agreed with Guintivano et al.'s estimates for prefrontal cortex. Our CHARM cell proportion estimates are on average higher than those obtained using 450K arrays, as the CHARM data were sampled using 2 mm dermal biopsy punches to minimize white matter contamination. The mathematical reconstruction of the methylation signal was also done for the unsorted samples (Additional File 1, Figure S4c).

Given that sorted data on the 450K array are only available for one brain region, we cannot demonstrate our improved ability to detect true brain-region differences in cell-type specific methylation on this platform. However, to show our ability to reduce false-positive signal, we constructed an artificial comparison by grouping the mixture samples with the highest and lowest neuronal fractions and applied models M1 and M2 to look for differences between these two groups. Any such differences are clearly due only to cell-fraction variation, and model M2 reduces the number of false-positive signals (Additional File 1, Figure S4d), as we saw for our CHARM data (Figure 2a). These results indicate that our methods apply well to data from the 450K array.

Generative model of methylation signal

To illustrate our model, we consider the case of estimating differences in methylation between DLPFC (D) and HF (H). We assume these brain tissues are composed of two cell types, NeuN+ (+) and NeuN- (-). For a fixed genomic position, we let μj,k be the methylation level in region j, j ∈ {H, D} and cell type k, k ∈ {+, -}. Scientifically, we are interested in identifying genomic locations where μH,k - μD,k ≠ 0, that is, where NeuN+ or NeuN- have different methylation levels in the two brain regions.

then our estimated coefficients have interpretations equivalent to the generative model in Equation 1. Specifically, we can test the hypothesis of no difference in NeuN+ methylation between D and H (μ_H,+ - μ_D,+ = 0) by testing the hypothesis that β₂ = 0, and the hypothesis of no difference in NeuN- methylation between D and H (μ_H,- - μ_D,- = 0) by testing the hypothesis that β₃ = 0.

From the equations above, we can see that estimating the fraction of cells of each type, π_i, allows us to explicitly find locations with brain-region differences specific to NeuN+ or NeuN- cells.

Naïve models are biased

This means that where + and - have different methylation levels (δ ≠ 0), a difference in the fractions of + and - cells in the different brain regions can lead to false-positive signals of brain region differences in methylation.

With $\bar{\pi }$ and ${\bar{\pi }}_{{}^H}$ the average of the ${\pi }_i$ 's in all samples and H samples, respectively, and ${\bar{\pi }}^2$ and ${\bar{\pi }}_H^2$ the average of the ${\pi }_i^2$ 's in all samples and H samples, respectively. Note that the bias is directly proportional to the difference between NeuN+ and NeuN- fractions, demonstrating that this approach is incapable of deconvolving these quantities of interest.

Estimation of mixture proportions

If we assume ${\mu }_{D,+}$ and ${\mu }_{D,-}$ are known, ${\pi }_i$ is the only unknown in this equation, so it can be estimated. Note that we do need to constrain our estimate of ${\pi }_i$ to be between 0 and 1. Also, the means μ are not known, so we collected data to allow us to estimate these means, by measuring methylation in pure cell sorted + or - fractions from each brain region of interest. Given that these methylation measurements have uncertainty, we want to reduce the uncertainty in our estimate of ${\pi }_i$ by using many informative genomic regions. We first select a set of genomic regions where + and - methylation differs. We then find the optimal value of ${\pi }_i$ to explain the observed methylation for sample i in these locations, as a function of our estimated means and ${\pi }_i$, subject to the constraint that ${\pi }_i$ is between 0 and 1. This procedure closely follows that presented by Houseman et al. [19].

Selection of the genomic locations can be based on a variety of factors, such as the range of observed methylation at these locations, the variance of the methylation estimates, and the length of the region of differential methylation. For our estimation, we chose the 500 genomic regions which were the strongest + vs. - DMR candidates in the brain region of interest in relation to the amount of methylation difference and the length of the region showing the methylation difference. We found that our results were quite robust to the number of regions selected, with 500 performing well.

To investigate whether it is absolutely necessary to have sorted data from a given brain region to estimate cell proportions in unsorted data from that region, we identified a set of 'universal' genomic regions. These universal regions had different NeuN+ and NeuN- methylation signals within a brain region, but showed consistent NeuN+ and NeuN- methylation levels across the three brain regions for which we had data (DLPFC, HF, and STG). Many of these + vs. - DMR candidates had consistent NeuN+ and NeuN- levels across brain regions, with 14% to 17% of the probes in the + vs. - DMRs belonging to genomic regions of consistent signal. We estimated the means μ in these regions of consistent signal using sorted data from DLPFC alone, and then performed cell fraction estimation in the unsorted samples from DLPFC, HF, and STG using these mean values. Since we do not know the true cell fractions in these unsorted samples, we used the estimates we had obtained for each brain region using the region-specific DMRs and mean values, as described above, as our gold standard.

All analysis was implemented in R (R Core Team, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing: Vienna, Austria, 2012; [33]). The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO series accession number GSE48610.

Effect of inaccurate mixture estimates

where η_i is between 0 and 1, and the third line follows from the fact that γ_i must be between $-{\pi }_{\mathsf{\text{i}}}$ and $1-{\pi }_{\mathsf{\text{i}}}$ to ensure that ${\pi }_{\mathsf{\text{i}}}^{*}$ is between 0 and 1. We can see that the coefficient of $X_{\mathsf{\text{i}}}\left(1-{\pi }_{\mathsf{\text{i}}}\right)$ is no longer measuring just the quantity we are interested in (the difference between NeuN+ methylation in regions H and D), but it also has an additional factor related to the size of the estimation error, and similarly for the coefficient of $X_{\mathsf{\text{i}}}{\pi }_{\mathsf{\text{i}}}$.

CHARM DNA methylation analysis

Genomic DNA was isolated from brains using the Masterpure kit from Epicentre, according to the manufacturer's protocol. For genome-wide DNA methylation assessment, 1 ug of genomic DNA from each sample was digested, fractionated, labeled, and hybridized to a CHARM array as described [34, 35] using a custom Nimblegen 2.1 million feature array assaying 5,114,655 CpG sites. We used the Bioconductor package 'charm' for sample preprocessing along with the package 'bumphunter' for DMR identification and permutation computation.

Human postmortem brain samples

Fluorescence-activated cell sorting was performed on frozen postmortem dorsolateral prefrontal cortex (n = 4), and hippocampal formation (n = 4) and superior temporal gyrus (n = 3) from individuals not affected with neurological or psychiatric disease. To validate the statistical model, we used nine additional healthy samples from the dorsolateral prefrontal cortex. These samples underwent nuclei extraction and sorting as described below. The model was applied to additional unsorted control samples (19 samples from DLPFC, 13 samples from HF, 31 samples from STG) to deconvolve NeuN+ and NeuN- methylation signatures. All samples were obtained from the bank of the Center for Neurodegenerative Disease Research (CNDR) in the Department of Pathology and Laboratory Medicine at the University of Pennsylvania (directed by Dr John Q Trojanowski, see Additional File 1, Tables S2-4 for demographic information).

Nuclei extraction, NeuN labeling, and sorting

Total nuclei were extracted via sucrose gradient centrifugation as previously described [25]. A total of 250 mg of frozen tissue per sample was homogenized in 5 mL of lysis buffer (0.32M sucrose, 10 mM Tris pH 8.0, 5 mM CaCl₂, 3 mM Mg acetate, 1 mM DTT, 0.1 mM EDTA, 0.1% Triton X-100) by douncing 50 times in a 40-mL dounce homogenizer. Lysate was transferred to a 15 mL ultracentrifugation tube and 9 mL of sucrose solution (1.8 M sucrose, 10 mM Tris pH 8.0, 3 mM Mg acetate, 1 mM DTT) was pipetted to the bottom of the tube. The solution was then centrifuged at 27,000 rpm for 2.5 h at 4C (Beckman, L8-80 M; SW28.1 rotor). After centrifugation, the supernatant was removed by aspiration and the nuclei pellet was resuspended in 500 uL of PBS.

The nuclei were incubated in a staining mix (0.71% normal goat serum, 0.036% BSA, 1:1200 anti-NeuN NeuN (Millipore, MAB377), 1:1400 Alexa647 goat anti-mouse secondary antibody (Invitrogen, 21236) for 45 min by rotating in the dark at 4°C. Unstained nuclei and nuclei stained with only secondary antibody served as negative controls. The fluorescent nuclei were run through a FACS machine with proper gate settings. A small portion of the NeuN⁺ and NeuN^- nuclei were re-run on the FACS machine to validate the purity. Immunonegative (NeuN^-) nuclei were collected in parallel. To pellet the sorted nuclei, 2 mL of sucrose solution, 50 uL of 1 M CaCl₂, and 30 uL of Mg acetate were added to 10 mL of nuclei in PBS, incubated on ice for 15 min, then centrifuged at 3,000 rpm for 20 min. The nuclei pellet was resuspended in 10 mM Tris (pH 7.5), 4 mM MgCl₂, and 1 mM CaCl₂. Fluorescent images were taken on a Zeiss Axio Observer. Z1 microscope with a Plan-Apochromat 100x/1.40 oil-immersion objective lens. Images were generated using an Axiocam MR3 microscope camera and Axiovision software (AxioVs40, version 4.8.2.0, Carl Zeiss, Inc). Images were processed using ImageJ.