Accounting for cell lineage and sex effects in the identification of cell-specific DNA methylation using a Bayesian model selection algorithm

Cell- and sex-specific differences in DNA methylation are major sources of epigenetic variation in whole blood. Heterogeneity attributable to cell type has motivated the identification of cell-specific methylation at the CpG level, however statistical methods for this purpose have been limited to pairwise comparisons between cell types or between the cell type of interest and whole blood. We developed a Bayesian model selection algorithm for the identification of cell-specific methylation profiles that incorporates knowledge of shared cell lineage and allows for the identification of differential methylation profiles in one or more cell types simultaneously. Under the proposed methodology, sex-specific differences in methylation by cell type are also assessed. Using publicly available, cell-sorted methylation data, we show that 51.3% of female CpG markers and 61.4% of male CpG markers identified were associated with differential methylation in more than one cell type. The impact of cell lineage on differential methylation was also highlighted. An evaluation of sex-specific differences revealed differences in CD56+NK methylation, within both single and multi- cell dependent methylation patterns. Our findings demonstrate the need to account for cell lineage in studies of differential methylation and associated sex effects.


Introduction
DNA methylation is a widely studied epigenetic modification that plays an essential role in the regulation of gene expression [1], cell differentiation [2] and the maintenance of chromatin structure [3]. Advances in

Materials and methods Data
Publicly available Illumina 450K methylation data were obtained from six healthy males subjects [13]. The data contained methylation β-values on seven isolated cell populations: CD19B+ cells, CD14+ Monocytes, CD4 + T cells, CD8 + T cells , CD56 + Natural Killers (NK), Neutrophils (Neu) and Eosinophils (Eos). Samples of the same cell types excluding Eosinophils were also obtained for five healthy female subjects [18]. In the absence of Eosinophils, The methodology described in this paper was therefore only applied to the six cell types common to both datasets. For both sexes, 445,603 CpG probes were available for analysis.
Consistent with previous findings [13], differences in DNA methylation levels across CpG sites were primarily driven by cell type differences, as opposed to being an artefact of between-subject variation (Fig 1).
A hierarchical clustering of cell-sorted samples revealed that differences in methylation levels were closely aligned with expected cell lineage in whole blood and were sex-independent. Additional cell sorted 450K methylation data on female and male subjects [18], aged 32-50 years, were obtained from public databases, for validation of model results. Female data consisted of cell-sorted samples on six healthy subjects downloaded from ArrayExpress (Accession number: E-ERAD-179) for CD19B + B cells, CD14 + Monocytes, CD4 + T and CD8 + T cells. Methylation data on six healthy male subjects from the same study were downloaded from Gene Expression Ominibus (Accession number: GSE71245). Following filtering and normalisation, 436,067 and 445,307 CpG sites were available for validation for females and males respectively, based on their correspondence with available sites in the training data. Table 1: Description of cell-specific models based on cell lineage from Fig 1. Each candidate model corresponds to a partition of cell types into non-overlapping groups. Cell types within each set of parentheses, {}, belong to the same partition and were assumed to have the same level of methylation. For each cell-specific model excluding 'All', the partition annotated by an asterisk ( * ) denotes the reference partition.

Model Formulation
The methodology presented in this paper was motivated by the identification of cell-specific methylation at the individual CpG level. A cell-specific CpG was defined as any CpG probe where the expected methylation level was different for one (single cell-specific) or more (multi-cell-specific) cell types, relative to others observed. The true cell-specific methylation pattern at each CpG probe was assumed unknown and treated as a model selection problem, with key details outlined in this Section. Information about cell lineage inferred from available samples was used to define candidate models for selection at each CpG. Using the results of hierarchical clustering (Fig 1), eleven candidate models were defined (Table 1) where each model represented a dendrogram split or node. Each candidate model therefore proposed a different cell-specific methylation pattern, characterised by a unique partition of cell types into non-overlapping groups. For a given candidate model, cell types assigned to the same group were assumed to share the same level of methylation. Two additional models, corresponding to the null and saturated cases, were also considered.
where N J (a, A) defines a J−dimensional Normal distribution with mean vector a and variance-covariance matrix A. The unknown variance was assumed to be common across all cell types, denoted by σ 2(m) ks I J×J where I is the identity matrix. This variance-covariance structure was chosen for identifiability reasons given sample sizes available.
A Bayesian approach to model selection was adopted, which allowed for probabilistic statements to be made about the relative fit of each candidate model. Under this approach, the posterior probability of each model conditional on the observed data was calculated for all candidate models. The posterior probability of model m compared to other candidates m ′ was calculated as, P r (Model m|y ks ) = p (y ks |Model m) P r (Model m) where the sum of probabilities over all candidate models was equal to 1. The term p (y ks |Model m) was obtained by integrating over all unknown parameters from the likelihood and prior distributions. Prior probabilities of each model, P r (Model m), were assumed to be equal to reflect a lack of model preference a priori.
Prior distributions for remaining parameters were selected such that Eq 1 could be derived analytically.
A g-prior distribution [19] was adopted for each µ The g-prior is a popular choice in linear model selection settings, as it allows the experimenter to introduce information on the scale of X m . Expected methylation levels were centered around the overall methylation level, b 0 , with variance proportional the standard error of each partition. The scaling factor g s > 0 is interpreted as a relative weighting of the prior versus the observed data. Here, g s was assumed common over all CpG probes and estimated by maximising the marginal likelihood averaged over candidate models. In this paper, b 0 was set to the global methylation mean, averaged over CpGs and cell types. A conjugate prior distribution for σ 2(m) ks of form p σ 2(m) ks ∝ 1/σ 2(m) ks completed model specification [19]. The expectation-maximisation (EM) algorithm was used to obtain a global Empirical Bayes estimate for g s [20]. This computational approach provided significant computational benefits over sampling-based approaches, namely Markov chain Monte Carlo (MCMC), and was possible given closed-form solutions for p (y ks |Model m) conditional on g s . In addition, the desired posterior model probabilities from Eq 1 were available upon convergence of the EM algorithm. The proposed methodology was implemented in R, with code available as Supporting Information (File S1).

Making the most of the Bayesian approach for model based inference
The accommodation of model and parameter uncertainty under the Bayesian approach formed the basis of subsequent inference, namely the identification of cell-specific CpGs, the estimation of differential methylation by cell type and the assessment of sex-specific differences by cell type. Brief details of each inference are provided in this Section.

CpG marker identification
CpG probes or 'markers' associated with each cell-specific pattern were identified by comparing posterior model probabilities at each CpG probe. For each candidate model and sex, markers were identified by reviewing the set of K model probabilities: P r (Model m|y 1s ) , . . . , P r (Model m|y Ks ). A 5% Bayes' False Discovery Rate (FDR) [21] was applied to each set of probabilities to control the expected number of false discoveries. Common CpG markers were defined as CpGs that were identified in both female and male samples, for the same candidate model. The compilation of common marker panels allowed for the assessment of sex-specific differences by cell type, within each candidate model.

Estimation of differential methylation
Posterior inference about each marker panel focused on the estimation of differential methylation by cell type relative to its corresponding reference partition (Table 1). Under model m, the posterior distribution of differential methylation for cell-specific partition p and sex s is, where c p is the number of cell types in partition p and µ The posterior distribution in Eq 2 was also used to validate selected CpG markers. Validation of common CpG markers was limited to CD19 + B, CD4 + T, CD8 + T and Pan T models, in light of validation samples available. Using all validation samples for each sex, the average β−value for each cell type and CpG was computed. The difference between each average β−value and reference partition was then compared with the corresponding 95% CI from Eq 2 for the appropriate candidate model. Concordance rates with respect to the predicted direction of methylation (hypomethylated, hypermethylated) for validation samples were also calculated. This joint approach was motivated by the limitation that validation based on coverage of 95% CIs relied on a posterior estimate of σ 2(m) ks . When this estimate is small and/or underestimated, proportions of validated markers based on CI coverage only were likely to be low, even if differences in methylation between training and validation samples were small. Finally, it is noted that not all cell types observed in the training data were available in the validation samples. Given the properties of the multivariate Normal distribution, this discrepancy was addressed by using the appropriate marginal distributions for the available cell types.

Evaluation of sex effects
A similar expression to Eq 2 was derived to evaluate sex-specific methylation differences within common CpG marker panels. In this case, focus was on the comparison of female and male methylation estimates for differentially methylated cell types, as defined by the corresponding candidate model. The posterior distribution of this difference across model partitions was Multivariate Normal, To assess evidence for sex effects, the posterior probability that the difference in methylation between sexes using Eq 3 was at least 0.10 was calculated for each cellspecific partition. This calculation again relied on estimates of the residual variance, σ 2(m) ks for each sex which were set to their posterior mean estimate. A sex-specific difference for a given partition was declared if the posterior probability of a difference greater than 0.10 exceeded 0.95.

CpG markers to genes: Assessment of SNP effects, genomic features and pathway enrichment
To provide additional evidence to support our method of marker identification, a pathways enrichment analysis was performed to explore the underlying biology of gene sets derived from common CpG marker panels. Common CpG markers were mapped to genes using Illumina Human Methylation 450k annotation data available from Bioconductor [22]. SNP information was also collated to infer the percentage of SNP associated markers with associations limited to SNPs located directly on the CpG loci. Using KEGG functional analysis in WebGestalt [23], a hypergeometric test was applied to each marker panel. A 5% FDR [24] was applied to resulting p-values to identify significant pathway enrichment for derived gene lists.

Cell lineage impacts the identification of differentially methylated CpGs
Across all cell-specific models, 83,449 and 97,747 CpG markers were identified for females and males, respectively ( Table 2). Among female samples, 42,834 markers (51.3%) were associated with differential methylation in multiple cell types, which included a three-fold increase in Pan T markers compared with males. A larger proportion of multi-cell dependent markers among males was observed (64.8%), due to larger numbers identified under Lymphocyte-II and Myeloid models. Among single cell dependent models, CD19 + B was the most frequently observed marker type for both sexes, with 25,611 in females and 18,271 in males.
A total of 42,452 CpGs were identified under the same cell-specific methylation pattern for both sexes, corresponding to 9.5% of the observed methylome. Within this subset, 23,551 (55.5%) were defined by differential methylation in more than one cell type. Over all candidate models, smaller frequencies of common markers were associated with T-cell dependent markers (CD4 + , CD8 + , Pan T).

Differential methylation is affected by cell lineage among common CpG markers
Common CpG markers associated with differential methylation in Myeloid cell types (CD14 + Mono, CD16 + Neu) were consistently hypomethylated across all relevant marker panels (Table 3). Among Lymphocytes, CD8 + T markers were the least likely to be hypomethylated (35.99%). Smaller proportions of hypomethylation among Lymphocyte I and II panels were indicative of greater methylation among Lymphocyte versus Myeloid cell subtypes.
The impact of cell lineage on differential methylation was greatest among marker panels related to Lymphocytes (Fig 2). Lower levels of differential methylation (<0.10) were concentrated within single cell   The distribution of differential methylation by sex among common Pan T markers revealed considerable variation for both hypermethylated and hypomethylated states compared to CD4 + T and CD8 + T panels (Fig 3). Differential methylation levels among Pan T markers tended to be greater for CD4 + T cells; 25 markers in this subset showed strong evidence of differential methylation greater than 0.5 for both sexes ( Fig  S1).

Validation of common CpG markers for B-and T-lymphocytes
The validation of common markers with respect to mean differences was generally higher in males than females, across all immune cell subtypes (Table 4). Whilst validation based on coverage of credible intervals was low, differences between training and validation estimates were relatively small across all markers tested; for CD4 + T markers, approximately 80% of absolute differences were less than 10% (Fig S2). Furthermore, comparisons with respect to inferred methylation state showed very high concordance rates for all marker panels and these findings were consistent for both sexes.

Genomic feature distributions and enrichment analysis
The effects of cell lineage were evident in the comparison of genomic features, with higher proportions of markers residing in Transcription Start Sites (TSS) for Lymphocyte cell subtypes compared with Myeloids (Table 5). Common markers were concentrated in the gene body and intergenic regions, with 54.73% of CD56 + NK markers located in the gene body. Among Lymphocyte cell subtypes, TSS proportions were highest for T cells, with 18.28% and 16.48% for CD4 + T and CD8 + T cells, respectively. Moderate proportions of markers were directly associated with SNPs and these levels were maintained between sex-specific and common marker panels, averaging 30.7% in common markers (Table S1). Higher SNP proportions were observed in marker panels related to the lineage of Myeloid cells, except for CD56 + NK markers of which 36.02% were SNP associated. For CpG probes not assigned to any common marker panel, a similar degree of SNP association was observed (26.61%). Significant pathways enrichment for common markers were associated with immune cell subtypes (Table   6) all including biologically relevant pathways.
Assessment of common marker profiles shows sex effects in CD56 + NK, CD16 + Neutrophils The majority of common CpG markers did not exhibit sex-specific differences in methylation profile (Table 7).  markers exhibited sex effects, with the majority corresponding to autosomal CpGs. No sex-specific differences were identified among common CD8 + T markers.
Sex-specific differences identified among single cell-specific markers were uniquely mapped to 215 genes (Table S2). The majority of identified cases were concentrated in the CD56 + NK common marker panel.
Four CD4 + T markers were associated with greater methylation observed in females and were mapped to a single gene, CD40LG, located on the X chromosome. Annotation information for the full list of sex-specific  Table 7: Summary of sex-specific differences by common CpG marker panel. A difference between male and female methylation estimates ≥0.10 was the outcome of interest. A marker was declared sex-specific if the posterior probability for this outcome exceeded 0.95, for one or more model-based partitions. For each marker panel, total numbers of sex-specific markers and autosomal sex-specific markers are given.
Non  (File S2). sex effects in CD56 + NK methylation was prominent in Lymphocyte-I and Lymphocyte-II common marker panels, in addition to the CD56 + NK panel (Fig 5). Within the Lymphocyte-I panel, 1698 common CpG markers were identified as sex-specific, of which 1627 showed differences in CD56 + NK methylation only. Similarly, 442 common Lymphocyte-II markers showed differences in CD56 + NK methylation between sexes. The tendency for CD56 + NK methylation to be higher in males was common across all three panels. 684 common Myeloid markers were associated with differences in CD14 + Monocytes or CD16 + Neutrophils only. Differences with respect to CD14 + Monocytes tended to be higher in females, compared with higher male methylation in CD16 + Neutrophils (Fig 6).

Discussion
This paper has proposed new statistical methodology for the discovery of cell-specific methylation profiles in whole blood, applying principles of Bayesian model selection. The characterisation of CpGs by differential methylation in one or more cell types builds significantly upon existing work that has been restricted to univariate analyses of differential methlyation by cell type. Sex-specific differences in both the prevalence of cell-specific profiles and methylation signal for select cell types were also demonstrated. For immune cell subtypes, validation of common markers using external cell-sorted samples produced favourable results. Enrichment analyses of common marker panels provided additional support for the proposed methodology, where it was demonstrated that the detection of cell-specific methylation at the individual CpG level was also biologically meaningful.
The incorporation of knowledge about hematopoietic lineage into model specification represents a semisupervised approach that reflects relevant cell biology. Consistent with other model selection strategies, our approach assumes that the defined set of candidate models is exhaustive and true patterns beyond this set are not of primary interest. In the event that the true methylation pattern is not commensurate with any of the candidate models specified, two outcomes are likely. In the first instance, corresponding probes are It is well documented that there are sex-specific differences in the proportions of circulating white blood cells [25,26]. The application of the proposed methodology to female and male samples has highlighted the importance of accounting for sex effects in DNA methylation analyses. Greater numbers of CD19 + B and T-cell dependent markers in females are consistent with previous findings and are possibly indicative of higher levels of cell activation [27]. The association between sex-specific differences in select CD4 + T markers and the CD40LG gene have also been identified previously. Previous studies have pointed to alelle specific methylation for this gene [28,29] where CD4 + T hypomethylation is observed in healthy males compared with healthy women who carry one methylated and one hypomethylated alelle. One of our major findings was large differences of methylation between males and females in markers defined as CD56 + NK specific. This is interesting when considered alongside the observation that males show an increase in circulatory NK cells compared to females [27], which adds further support for the accuracy of the approach. Additionally there is some evidence of sex-specific methylation differences in CD + 56 NK, as well as CD + 8 T-cells [30]. Under the proposed approach, we have provided a potential solution to accurately account for potential bias introduced by sex effects at the marker level.
The presence of cell-specific methylation CpG markers highlights the need to account for cellular composition prior to conducting Epigenome Wide Association Studies (EWAS), in whole blood. Methods for this purpose have been developed [31,32] based on the assembly of methylation 'signatures' from cell-sorted data which are then projected onto heterogeneous samples to predict cell type proportions. A comparison of common markers with the top 500 CpG probes identified by the cell mixture methodology of [31] revealed 75% concordance between panels (data not shown), of which the majority were associated with Lymphocyte-I/II and Myeloid maker panels. The inclusion of other marker panels in these algorithms may lead to further improvement in cell mixture estimation, in particular for immune cell subtypes that may be present in low proportions. Furthermore, the performance of these algorithms rely on consistent cell-type effects across cohorts [33]. Given the sex-specific methylation differences we have identified in this study, failure to account for sex effects may also impact upon the quality of cell mixture estimation and should therefore be given due consideration.
It is common practice in array-based methylation studies to exclude CpG sites which contain SNPs both within the probe and on the CpG site. While this is a valid approach to filtering before analysis, it will often lead to dramatic reduction of overall data. As a result, it is likely that sites of potential interest may be lost before any association can be made. By mapping hg38 annotated SNPs to all 450K CpG loci, we were able to ascertain the overall proportion of cell marker sites which have a SNP present; on average, across the common set of markers, this was approximately 30.7% of markers. In light of these results, we suggest that deconvolution studies and methods should account for SNP events at cell marker sites, noting the proportion that are present. For the filtering stage, we recommend that the overall rarity of the SNP variant be taken into account, for example, retaining CpGs which also have a 'rare' (MAF < 0.01) variant mapping. This approach is likely to be beneficial to the overall study design and outcome.

Supporting information
Supplementary files S1 and S2 are available from the GitHub respository https://github.com/nicolemwhite/BayesMS. S1 File R code for Bayesian model selection algorithm. Core R functions required to prepare data and compute posterior model probabilities via the EM algorithm, across all listed candidate models.
S2 File List of common CpG markers associated with a sex-specific difference of ≥0.10. Illumina Human Methylation 450k annotation data are also included for each CpG marker.  Table S2: Distribution of sex-specific common markers over chromosomes, for single cell-dependent makers only. A sex-specific marker was declared if the posterior probability from Eq 3 was greater than 0.95 for at least one differentially methylated cell type.