BACKGROUND
dna methylation  is one of the most important epigenetic mechanisms regulating gene expression, and aberrant dnam has been implicated in the initiation and progression of human cancers  <cit> .
dnam changes have also been observed in normal tissue as a function of age  <cit> , and age-associated dnam markers have been proposed as early detection or cancer risk markers  <cit> .
proper statistical analysis of genome-wide dna methylation profiles is therefore critically important for the discovery of novel dnam based biomarkers.
however, the nature of dna methylation data presents novel statistical challenges and it is therefore unclear if popular statistical methods used in the gene expression community can be translated to the dnam context  <cit> .
the illumina infinium humanmethylation <dig> beadchip assay is a relatively recent high-throughput technology  <cit>  that allows over  <dig>  cpgs to be assayed.
while a growing number of infinium 27k data sets have been deposited in the public domain  <cit> , relatively few studies have compared statistical analysis methods for this platform.
in fact, most statistical reports on infinium 27k dnam data have focused on unsupervised clustering and normalisation methods  <cit> , but as yet no study has performed a comprehensive comparison of feature selection and classification methods in this type of data.
this is surprising given that feature selection and classification methods have been extensively explored in the context of gene expression data, see e.g.
<cit> .
moreover, feature selection can be of critical importance, as demonstrated by gene expression studies, where for instance use of higher order statistics has helped identify important novel cancer subtypes  <cit> .
given that the high density illumina infinium 450k methylation array is now starting to be used  <cit>  and that this array offers the coverage and scalability for epigenome wide association studies   <cit> , it has become a critical and urgent question to determine how best to perform feature selection on these beadarrays.
the illumina infinium assay utilizes a pair of probes for each cpg site, one probe for the methylated allele and the other for the unmethylated version  <cit> .
the methylation level is then estimated, based on the measured intensities of this pair of probes.
two quantification methods have been proposed to estimate the methylation level:  β-values and  m-values.
while the β-value measures the percentage of methylation at the given cpg site and is the manufacturer's recommended measure to use  <cit> , the m-value has recently been proposed as an alternative measure  <cit> .
the m-value is defined by the log <dig> ratio of the intensities of the methylated to unmethylated probe and has also been used in the context of other methylation array technologies  <cit> .
the m-value is also an analogue of the quantity which has been widely used in expression microarray analysis, although there are two important differences: in the infinium dnam  assay, both probes are  always measured in the same colour channel, thus dye bias does not need to be adjusted for, and  the two probes are measured in the same sample.
a recent study based on a titration experiment compared β and m-values and concluded that m-values, owing to their more homoscedastic nature , were a better measure to use  <cit> .
however, this study was limited to one data set and feature selection was only investigated using fold-changes, while statistics and p-values were not considered.
the purpose of our study is therefore two-fold:  to provide a comprehensive comparison of the performance of β and m-value measures in the context of popular feature selection tools that use actual statistics to rank features , and  to compare the performance of some of the most popular feature selection and classification methods on dnam data from different tissues and correlating with different phenotypes, so as to gain objective insights into how best to perform feature selection and classification in this novel context.
to address these goals, we make use of  <dig> independent sets of human dna methylation data, all generated with the 27k illumina infinium platform, encompassing over  <dig> samples and representing over  <dig> million data points.
materials and methods
data sets
the  <dig> human dna methylation data sets are summarised in table  <dig>  all the data were generated using the illumina infinium humanmethylation <dig> beadchip assay that enables the direct measurement of methylation at over  <dig>  individual cytosines at cpg loci located primarily in the promoter regions of  <dig>  unique genes.
all data sets, except the tcga lung cancer set, followed the same quality control and normalisation strategy.
briefly, non background-corrected data was used and only cpgs with intensity detection p-values less than  <dig>   were selected.
inter-array normalisation correcting for beadchip effects and variations in bisulfite conversion efficiency were performed using a linear model framework with explicit adjustment for these factors, but only if these factors were significantly correlated with pca/svd components.
the  <dig> independent 27k dnam data sets used in this paper: ukops - uk ovarian cancer population study: whole blood samples from women with ovarian cancer and age-matched healthy controls; endom - normal endometrium and endometrial cancer; cervx - normal cervix and cervical cancer; ovc - ovarian cancer tissue; t1d - type  <dig> diabetes: whole blood samples from type- <dig> diabetics; bc - normal breast tissue and breast cancer; and lc - normal lung tissue and lung cancer.
number of cpg probes passing quality control, total number of samples, tissue type, number of normal/cancer samples, age-range and reference to data are given.
na: not available.
definition of β-value and m-value
in the illumina infinium human methylation 27k beadchip assay, bisulphite  converted dna is amplified, fragmented and hybridised to the beadchip arrays, with each chip accommodating  <dig> samples as designated by sentrix positions a-l. each interrogated locus is represented by specific oligomers linked to two bead types, with one representing the sequence for methylated dna  and the other for unmethylated dna .
the methylation status can be measured either by β-values or m-values.
for each specific cpg site, the β-value is calculated from the intensity of the methy and unmethy alleles, as the ratio of the fluorescent signals:
β=maxmethy,0/maxmethy,0+maxunmethy,0+ <dig>
the β-value is a continuous variable between  <dig>  and  <dig>  representing the ratio of the methylated allele to the combined locus intensity.
in contrast, the dna methylation m-value is calculated as the log <dig> ratio of the intensities of the methylated probe to the unmethylated probe:
m=log2maxmethy,0+1/maxunmethy,0+ <dig>
the m-value is therefore a continuous variable which can in principle take on any value on the real line.
since most interrogated cpg sites have intensities  larger than  <dig> , the offset values  have only a small effect on both measurements.
hence, the relationship between the β and m-values is approximately logistic:
β=2∧m/2∧m+1;m=log2β/1-β.
while the β-value has a more intuitive biological interpretation it suffers from severe heteroscedasticity   <cit> .
in contrast, the m-value is not directly interpretable in terms of a methylation percentage, but is more statistically valid for analysis of differential methylation levels owing to its more homoscedastic nature  <cit> .
evaluation of feature selection methods
we use a multiple training-test set strategy  <cit>  to compare the true positive detection rates of three different feature selection methods:  without filtering ,  filtering based on variability ,  filtering using supervised principal components .
our justification for focusing on vf and spca is that both methods have been very popular and effective in the context of gene expression studies  <cit> .
the performance measure we use to compare the different feature selection methods is the estimated positive predictive value .
that is, for each method we select features significantly associated with a phenotype of interest using a training set, and these features are then evaluated in the test-set.
those that are also significant in the test set are called true positives, so the ppv measures the probability that a feature declared to be positive in the training set is a true positive.
we point out that strictly speaking a true positive can only be a feature which has been shown to be positive using a gold-standard technology.
however, in the absence of such gold-standard data, the training/test set strategy advocated here provides a suitable evaluation framework, since a feature that cross-validates is more likely to represent a true positive than one that does not.
moreover, this strategy has already been successfully used in previous work, see e.g.
<cit> .
thus, for each of the data sets in this study,  <dig> independent training and test sets  were randomly generated ensuring similar demographics for the phenotype of interest in training and test set.
for each training set, a supervised analysis is used to evaluate the association between each cpg methylation profile  and the phenotype of interest.
in the case of a continuous phenotype  the association is carried out under a linear regression model framework, while for a binary phenotype  we use t-statistics.
test statistic, p-value, q-value  <cit> , absolute difference in means, and variance of each cpg are recorded.
variance-filtering  is a feature selection method that filters the cpgs based on their variability before the supervised analysis is performed  <cit> .
the steps are
i. in each training set, select the  <dig>  most variable cpgs;
ii.
rank the  <dig>  most variable cpgs according to the significance of their p-values  and retain the top ranked  <dig>  cpgs .
.
iii.
using stepwise evaluation sets , count the number of true positives: i.e.
how many cpgs have a statistic p-value <  <dig>  in the corresponding test set while exhibiting the same directional change  as that in the training set;
iv.
for each evaluation set, the positive predictive value  is defined as the percentage of cpgs that are true positives.
supervised pca  is a powerful semi-supervised approach for feature selection and classification, first proposed in  <cit> .
here we adapt it to the feature selection context.
the steps are:
i. in each training set, perform the supervised analysis and rank all cpgs according to the significance of their p-values;
ii.
select the top ranked  <dig>  cpgs;
iii.
run singular value decomposition  analysis on the top ranked  <dig>  cpgs and find the principal component most correlated to the phenotype of interest;
iv.
rank the  <dig> cpgs according to the absolute value of their corresponding coefficient in the selected principal component and further select the  <dig>  cpgs with the largest absolute weights;
v. use the stepwise evaluation sets and count, for each evaluation set size, the number of true positives in the test set using the same criteria as for the vf method;
vi.
for each evaluation set, compute the ppv as described for vf method.
in the without filtering  method we perform for each training set the supervised analysis on all ~ <dig>  cpgs , select the top  <dig>  and then perform the same stepwise evaluation set analysis to estimate ppvs from the p-values in the corresponding test set.
the result of applying each of these  <dig> methods is a matrix of ppv values with rows labelling the evaluation set size  and columns labelling the  <dig> different runs .
choice of evaluation set size
given that different studies may exhibit widely different operating characteristics , even in relation to the same phenotype of interest, we decided that a more objective comparison of different feature selection methods across studies would be facilitated by choosing evaluation set sizes in each study that result in broadly similar ppv values across studies.
hence, evaluation set sizes were chosen so as to minimise differences in ppv between studies, subject to reasonable constraints on the minimum and maximum set sizes.
comparison of classification algorithms
we use the same multiple training-test partition strategy  as used in our feature selection comparison, to evaluate different classification algorithms.
thus, within each study, training and test sets were always of the same size and were balanced for the phenotype of interest.
we focused on four different powerful classification algorithms, which have been popular in the gene expression field:  supervised pca   <cit> ,  the lasso algorithm  <cit> ,  the elastic net   <cit>  and  support vector machines   <cit> .
we used the implementations of these algorithms as provided in the r-packages superpc, glmnet and e <dig>  in the case of spca we considered classifiers built from  <dig> up to  <dig> principal components, and the number of features  was optimised in the training sets using internal 10-fold cross-validation.
lasso is a special case of the elastic net with the penalty parameter α =  <dig>  <cit> .
the elastic net itself was run with α =  <dig>  .
the additional penalty elnet parameter λ was estimated in the training sets using internal 10-fold cross-validations.
the svm was run with a radial kernel with parameter γ =  <dig>  and with error term ε =  <dig> .
these svm parameters were fixed and not optimised, since we wanted to compare the algorithms at the same level of computational complexity .
we focused on age as the phenotype of interest since it represents a challenging scenario of small effect sizes and yet it is also well established that age has a significant impact on dna methylation patterns  <cit> .
since age is a continuous variable, svm was run in regression mode, and thus for all methods the predictor is a continuous score.
to evaluate concordance of age with the predictor in the test sets we used the c-index .
comparison of dimensional reduction algorithms
we considered two popular dimensional reduction algorithms: singular value decomposition   <cit> , and non-negative matrix factorization   <cit> .
in this work we use the svd implementation as computed by the lapack/linpack routine available in r http://www.r-project.org.
the application of nmf to dna methylation is justified due to the positivity of beta-valued data.
to perform nmf, the "nmf" r package  <cit>  was used.
nmf was run using the "brunet" algorithm and initialised using non-negative double svd .
hence, using this nndsvd initialisation nmf yields the same solution under repeated runs and allows for a direct comparison to svd.
the evaluation of the  <dig> algorithms was performed by comparing the correlations of the inferred components with the phenotypes of interest , specifically we derived and compared the corresponding r <dig> values under a linear regression model.
in the case of nmf, we needed to specify the number of components to infer.
this number was estimated from the svd analysis by comparing the spectral eigenvalues to the corresponding ones obtained by randomly permuting elements in the data matrix.
this number was also the number of svd components used when evaluating correlations with the phenotype of interest.
RESULTS
the signal to noise and signal strength landscape of dna methylation studies
one of the aims of this study is to compare different feature selection and quantification methods in relation to important study-specific parameters.
one of these parameters is the typical effect size of the study.
in the case of a binary phenotype, the effect size of a cpg is defined, loosely speaking, as the ratio of the difference in means between the two phenotypes to a pooled standard deviation.
in other words, it can be thought of as a signal to noise ratio, where the noise captures both biological and technical variation.
another important parameter to take into consideration in any association analysis is the number of features with an effect size larger than some significance threshold.
thus, in this study we ask if the effect size of cpgs associated with a phenotype of interest  and their number  have an impact on the performance of the different feature selection methods and if this depends on the methylation measure  used.
in order to consider a wide range of different effect sizes and signal strengths, we considered two main phenotypes: cancer/normal status  and age; and two different tissue types: epithelial tissue and whole blood.
the corresponding 27k data sets used in this study are summarised in table  <dig>  we verified that these studies exhibited a wide range of different effect sizes and signal strengths, depending largely on tissue type .
for instance, we observed, as expected, that cancer diagnostic dnam markers in epithelial tissue have large snrs and signal strengths, in contrast to cancer diagnostic markers in blood, which, while numerous, were characterised by much smaller effect sizes .
this fits in with the expectation that cancer associated changes in whole blood reflect mostly changes in the underlying blood cell type composition  <cit> .
in stark contrast to diagnostic markers, we observed that age has a relatively much weaker impact on dnam patterns , and therefore estimated effect sizes were also smaller .
the signal strength of age-associated cpgs also varied depending on the tissue in which they were measured .
the signal strength  and the average effect size  of the top  <dig> ranked cpgs for diagnosis and age-settings and for epithelial and blood tissue.
beta-values were used to estimate effect sizes.
in the case of age, samples were divided into two groups by the median.
m-values significantly outperform β values but only in the small sample size limit
a previous study based on a titration experiment advocated the use of m-values over β values, since m-values exhibit less heteroscedasticity  <cit> .
however, this study only used fold-changes to rank features.
it is therefore unclear if results differ had the features been ranked using statistics.
thus, we asked if using β-values or m-values to quantify methylation has a differential impact on the statistics of differential methylation and on the positive predictive value .
first, we verified that in our data sets, β values were highly heteroscedastic, with very low variability at the extremes of the beta-scale, while m-values were approximately homoscedastic .
using t-statistics and q-values to estimate the false discovery rate  at different p-value thresholds and using all available samples , we observed however that the difference in fdrs between β and m-values within each study was very minor .
next, we compared the ppvs obtained using either β or m-values using a multiple 50% training 50% test set strategy where features were ranked according to their t-statistics .
because for a given methylation measure we observed substantial differences in fdr rates between studies, in some cases even between studies that looked at the same phenotype of interest and tissue type, we decided to choose evaluation set sizes specific to each study to ensure that ppv values were as similar across studies as possible .
in most studies there was no appreciable difference in the ppv between β and m- values and the absolute magnitude of the t-test statistics of the top ranked cpgs were also very similar .
to study the dependency on sample size we considered again a multiple training-test set strategy but this time using training and test sets where the number of samples of each phenotype was low .
since ordinary t-statistics can not be used in this small sample size limit, we used instead a popular empirical bayes approach to rank features according to a regularized t-statistic  <cit> .
in this small sample size limit, we observed that the ppv was significantly improved when using m-values, in some instances by at least 10% .
interestingly, using 50% training/test partitions, corresponding to the largest possible sample sizes, the ppv for m-values was still higher than for β-values, but only by at most 2% .
however, for larger sample sizes, ppv values were close to  <dig> and thus differences between m-values and β-values would naturally be smaller.
thus, the improved performance of m-values over β-values is mainly due to the use of regularized t-statistics, since we also observed that in the large sample size limit ordinary t-statistics performed equally well irrespective of m-value or β-value basis .
we also observed that in the large sample size limit, ordinary t-statistics performed similarly if not better than regularized t-statistics .
thus, given that β-values are more directly interpretable and that in this work we are mostly interested in the performance for relatively large sample sizes, we henceforth focus mainly on β-values.
variance filtering  can significantly reduce the ppv in the setting of small effect sizes
next, we evaluated the effect of filtering cpgs based on variance prior to performing the supervised analysis.
variance filtering has been a very popular feature selection tool in gene expression studies , and it has been shown to improve the detection rate of true positives  <cit> .
hence, we asked if the same result holds in the context of dna methylation data.
we observed that filtering based on variance improved the ppv in  <dig> of the  <dig> cancer/normal studies and in  <dig> of the  <dig> studies when age was the phenotype of interest .
in the age-setting, vf did worse than wf in two studies by approximately 5% in ppv values, while in the diagnostic setting the worse performance observed in two studies was far less marginal , indicating that in cases where effect sizes and signal strengths are small , filtering based on variance can be counterproductive.
thus, we posited that in the case of cancer diagnostic markers in whole blood , that vf would also lead to a decreased ppv performance if larger evaluation sets are used since the larger evaluation sets would include proportionally more false positives .
similarly, we would expect that for larger evaluation sets the wf method would allow proportionally more true positives to be captured.
confirming this, we observed that wf outperformed vf as the number of top ranked cpgs is increased .
this switch in performance between wf and vf as a function of evaluation set size is a consequence of the relatively small effect sizes combined with the relatively large signal strengths of diagnostic markers in whole blood.
spca outperforms wf and vf when effect sizes are small
supervised principal components has been shown to be a simple yet powerful algorithm for performing feature selection and classification in the gene expression context  <cit> .
we therefore decided to compare spca to the wf and vf methods in terms of feature selection.
interestingly, in the diagnostic setting spca did not always outperform vf or wf, although in all data sets it achieved very high  ppv values .
however, in the more challenging scenario of small effect sizes and signal strengths, as illustrated by using age as the phenotype of interest, we can see that spca achieved substantially better ppv values .
thus, principal components can identify true positives more reliably when signal strengths are small.
on the other hand, if signal strengths are large, a principal component may not have the plasticity to capture all the important true positives, and thus inevitably will include false positives among the less weighted features within the principal component .
we verified that in the diagnostic setting, spca ppv values were substantially higher and comparable to those of wf and vf when using only the top  <dig> cpgs in the selected pc , demonstrating that at more relaxed feature selection thresholds spca includes numerous false positives
signal correlation structure is improved using beta values
given that overall feature selection was optimal using spca, we next asked if performance would be similar had we used m-values.
remarkably, we observed that spca ppv values were generally higher when evaluated using beta values .
this was specially true in the diagnostic setting and was independent of tissue type, indicating that signal strength has a major impact on the inference of biologically relevant principal components.
interestingly, therefore, beta-values provide a basis in which the correlative structure of cpg methylation profiles associated with cancer status is better preserved.
in contrast, in the m-value basis, it appears that the correlation structure between biologically relevant cpgs is compromised leading to worse modelling of the biological variation.
the elastic net and svms outperform spca and lasso
next, we turned our attention to the classification task.
we considered a total of four different powerful classification algorithms, including spca, lasso regression, the elastic net  and support vector machines  .
we used the same training-test set partition strategy as with our feature selection analysis.
age was chosen as the phenotype of interest, since it is now well established that age has an impact on dnam patterns  <cit>  and because age-associated effects are of a relatively small magnitude, thus also providing a more challenging evaluation scenario for the classification algorithms  <cit> .
moreover, age-associated dnam markers may represent cancer risk markers  <cit> .
we considered the same studies as those used for the feature selection analysis in the age setting.
across all studies, we observed that either elnet or svm were superior to lasso and spca .
the fact that elnet outperformed lasso in all data sets extends previous results obtained on gene expression  <cit>  to the dna methylation context.
moreover, we observed that spca always performed optimally in the context of a single component , i.e.
adding more principal components to the classifier did not improve performance .
interestingly, spca- <dig> performed similarly to elnet and svm in the two smaller studies , suggesting that the rigidity imposed by principal components can be of an advantage in this setting.
unsupervised modelling of diagnostic and age effects
finally, we compared unsupervised algorithms in their ability to model diagnostic and age effects in dna methylation studies.
we considered two of the most popular unsupervised dimensional reduction algorithms from the gene expression field: singular value decomposition   <cit> , and non-negative matrix factorization   <cit> .
given that dna methylation beta-valued data are defined on the compact support  and are thus positively valued, it is justified to explore the application of nmf in this context.
to objectively compare nmf to svd, nmf was initialised using non-negative double svd, thus allowing us to determine if the positivity constraints of the nmf framework add biological value.
for both svd and nmf we computed the r <dig> correlations between the inferred components  and the phenotype of interest.
as a modelling criterion we compared the average r <dig> of the two best components, obtained from each algorithm.
focusing first on the diagnostic setting, we can clearly observe that nmf outperformed svd/pca in all  <dig> studies, including the study comparing blood samples from cancers and controls .
however, nmf was not conclusively better in the age-setting of much smaller effect sizes .
discussion
a previous study compared β to m-values in the context of a titration experiment using fold-change values to rank cpgs  <cit> .
the study concluded that m-values, owing to their reduced heteroscedasticity, provided a better basis in which to detect true positives.
this result is not unexpected because fold changes only consider the magnitude of differences in mean values, regardless of the underlying biological and technical variability.
therefore, using β-values and fold changes, true positive cpgs with low or high mean beta methylation values would not be highly ranked and would result in a significant proportion of false negatives.
as we have shown here, using many different data sets and using a multiple training test set strategy to estimate the ppv, t-statistics derived from m-values and β-values are comparable and lead to similar ppv values, but only in the scenario of relatively large sample sizes.
in contrast, when sample sizes are small and a regularized t-statistic must be used, m-values provided a significantly better basis for inference, sometimes by as much as 10% improvements in the ppv.
that regularized t-statistics perform relatively poorly in the β-value basis is not surprising, since the regularisation involves an empirical bayes approach in which a posterior variance is estimated from estimates specifying the prior, an approach which is known to be more sensitive to the precise distributional properties of the data and is therefore also more sensitive to deviations from normality.
therefore, our comprehensive study using statistics to rank features, indicates that the severe heteroscedasticity of β-values does not compromise the reliability of identifying true positives in the case of relatively large sample sizes, but that m-values offer significantly improved inference in the limit of small sample sizes.
a second important finding is that variance filtering does not always increase the ppv.
indeed, in scenarios where effect sizes are small but where there are numerous associated features, we have seen how variance filtering can compromise the detection of true positives when more relaxed significance thresholds are used.
this was illustrated in the context of cancer diagnostic methylation markers in whole blood, where the changes reflect underlying changes in blood cell type composition as demonstrated by us previously  <cit> .
in this setting, effect sizes are small but may involve a potentially large number of genes, indeed potentially all those genes whose expression differentiates blood cell types from each other and whose expression is under epigenetic regulation.
since very often there is a desire to perform gene set enrichment analysis  <cit>  on a set of top ranked features, thresholds must be chosen to ensure a relatively small false negative rate .
hence, in the scenario described here, variance filtering may be specially counterproductive since it leads to a larger fnr and hence could compromise the detection of enriched go-terms or pathways.
it follows that any decision to filter features based on variance must take into consideration the clinical and biological context of the study.
our third key result is the improved performance of spca in detecting true positives, specially in the more challenging scenario when effect sizes and signal strengths are small.
while spca has been shown to outperform many other feature selection methods in the context of gene expression data  <cit> , an analogous result was still lacking in the case of dna methylation.
our data therefore show that spca is a powerful feature selection algorithm, independently of the molecular profiles considered or their underlying statistical distributions.
the improved performance over vf also indicates that the correlative dnam patterns can be exploited to identify true positives more reliably.
interestingly, we observed that spca worked better in the context of β-values.
this may be surprising as one might be inclined to believe that components of maximal variation are more robustly identified in a basis that promotes variation, i.e.
using m-values.
however, it is also the case that m-values would aggravate the effect of outliers and therefore distort the true components of maximal biological variation.
therefore, illumina's recommendation to use β-values over m-values appears to be further justified by the fact that β-values, by virtue of being bounded, may provide a natural regularization and thus suppress the undesirable effects of potential technical outliers.
we also compared various powerful classification algorithms, including spca, elnet, lasso and svms.
our results clearly indicate that elnet and svm are superior classification methods in the case where effect sizes are small.
given that elnet also provides an automatic means of feature selection , it would appear to be a preferable choice over svm.
nevertheless, svm outperformed elnet in three studies, including the two smaller ones .
interestingly, we also observed that classification using spca performed optimally in the context of a single component of variation .
since multi-component spca  uses a training set to estimate an optimal linear combination of principal components that then makes up the predictor, it would appear that the estimated weights specifying the linear combination is not mirrored in the test set, hence why spca- <dig> is optimal.
it will be interesting to investigate if this result generalizes to other phenotypes where effect sizes may be small.
finally, we also compared two of the most popular unsupervised algorithms  in their ability to model differences between cancer and normal tissues and differences associated with age.
similar to the results obtained on gene expression, we observed that nmf also outperforms svd/pca in the context of dna methylation data.
indeed, the meta-genes inferred from the nmf were generally more highly correlated with cancer/normal status across all data sets considered.
hence, despite the increased computational complexity of nmf, this algorithm should be used whenever expected effect sizes are large.
on the other hand, if effect sizes are small  we did not find that the increased computational complexity of nmf offered any advantage over pca.
in the context of evaluating feature selection and classification methods, it is important to discuss the impact of potential batch effects.
indeed, as shown in many previous papers  <cit> , known and unknown batch effects can affect a substantial proportion of features in an experiment and lead to biased estimates of statistical significance.
in this work we used inter-array normalised data where beadchip effects and variations in bisulfite conversion efficiency were adjusted for using a linear model framework, as these factors were observed to account for significant components of variation in pca analysis .
we observed that not correcting for these factors may alter the absolute performance of any given method, yet the relative performance of the different methods was largely robust to whether inter-array normalisation was performed or not .
since the inter-array normalisation used only adjusts for known effects , we have not considered here the potential effects of unknown confounders.
to address this would require adapting surrogate variable analysis  approaches  <cit>  in the training/test set evaluation framework considered here.
this has only very recently made possible through an extension of sva, called frozen sva  <cit> .
thus, it will be interesting to further compare these feature selection methods in the context of the present dna methylation data sets, and thus further assess the potential impact of  batch effects on feature selection and classification.
while the insights obtained in this study have been derived from one particular array , it is likely that similar insights would apply to the recently released scaled up 450k infinium beadchip  <cit> .
we should point out however, that one of the key differences between the 27k and 450k platforms is that the probes on the 450k beadarray come in two different designs and are thus characterised by widely different statistical distributions  <cit> .
thus, the insights obtained here are likely to apply only to the data restricted to probes of one particular design.
in any case, the much higher density of the 450k array , means that feature selection in 450k data will be even more critical.
CONCLUSIONS
the optimal choice of methylation measure and feature selection method in the context of 27k dna methylation data is dependent on the sample size, expected effect sizes and signal strength of the specific study.
while spca is a powerful feature selection tool in dna methylation studies where effect sizes and signal strengths are small, the elastic net and svms generally outperform spca in the context of classification.
these insights are important for anyone embarking on large scale dna methylation profiling using 27k or 450k beadarrays.
abbreviations
dnam: dna methylation; pca: principal components analysis; svd: singular value decomposition; spca: supervised principal components; elnet: the elastic net; svm: support vector machines; nmf: non-negative matrix factorisation.
competing interests
the authors declare that they have no competing interests.
authors' contributions
jz and aet performed the statistical analyses.
mw contributed samples and obtained funding.
aet wrote the manuscript with contributions from jz.
all authors read and approved the final manuscript
supplementary material
additional file 1
typical scatterplots of mean methylation  vs standard deviation in methylation  using beta  and m values .
green line denotes a mean loess smoother.
a) the blood samples from the  <dig> healthy controls in the ukops study.
b) the  <dig> blood samples from the t1d study.
click here for file
additional file 2
false discovery rate  estimation table at different thresholds of top cancer diagnostic associated cpgs obtained from the supervised analysis  in the ukops, endom, and cervx data sets.
click here for file
additional file 3
the evaluation set size  is plotted against the mean ppv  for the different studies.
the mean ppv represents an average over the  <dig> training-test partitions and over the  <dig> feature selection methods .
a) diagnostic setting, b) age setting.
the selected evaluation set sizes in each study and for each phenotype of interest are marked in red.
note that in the case of age, differences in mean ppv between studies could not be minimized because of the additional constraint of a reasonable minimum set size.
click here for file
additional file 4
a-b) comparison of positive predictive values  obtained from beta  and m-values  using the without-filtering method and using  <dig> training-test set partitions .
a) diagnosis setting in ukops, endom, cervx, bc and lc with evaluation set sizes of top  <dig>   <dig>   <dig>   <dig>  and  <dig> diagnosis-associated-cpgs, respectively; and b) age setting in t1d, ovc, cervx, ukops, endom, and bc with the top  <dig>   <dig>   <dig>   <dig>   <dig>  and  <dig> age-associated-cpgs as evaluation sets, respectively.
c-d) comparison of absolute t-statistics obtained from beta and m-values.
c) diagnostic setting in ukops, endom, cervx, bc and lc; and d) age setting in t1d, ovc, cervx, ukops, endom, and bc.
all p-values shown are from a wilcoxon-rank sum test.
click here for file
additional file 5
mean ppv values between spca, wf and vf in the diagnostic setting but using only the top ranked  <dig> cpgs as evaluation sets.
click here for file
acknowledgements
this work was funded by the eve appeal and a grant from the uclh/ucl comprehensive biomedical research centre project  <dig> and undertaken at uclh/ucl, which received a proportion of its funding from the department of health nihr biomedical research centres funding scheme.
this work was also supported by the european network translational research in gynaecological oncology  of the european society of gynaecological oncology .
we thank allison jones, usha menon, helga salvesen and heidi fiegl for providing data.
aet is supported by a heller research fellowship.

 dna methylation is one of the most important epigenetic mechanisms regulating gene expression, and aberrant dnam has been implicated in the initiation and progression of human cancers. dnam changes have also been observed in normal tissue as a function of age, and age-associated dnam markers have been proposed as early detection or cancer risk markers. proper statistical analysis of genome-wide dna methylation profiles is therefore critically important for the discovery of novel dnam based biomarkers. however, the nature of dna methylation data presents novel statistical challenges and it is therefore unclear if popular statistical methods used in the gene expression community can be translated to the dnam context. the illumina infinium humanmethylation beadchip assay is a relatively recent high-throughput technology that allows over cpgs to be assayed. while a growing number of infinium 27k data sets have been deposited in the public domain, relatively few studies have compared statistical analysis methods for this platform. in fact, most statistical reports on infinium 27k dnam data have focused on unsupervised clustering and normalisation methods, but as yet no study has performed a comprehensive comparison of feature selection and classification methods in this type of data. this is surprising given that feature selection and classification methods have been extensively explored in the context of gene expression data, see e.g.. moreover, feature selection can be of critical importance, as demonstrated by gene expression studies, where for instance use of higher order statistics has helped identify important novel cancer subtypes. given that the high density illumina infinium 450k methylation array is now starting to be used and that this array offers the coverage and scalability for epigenome wide association studies, it has become a critical and urgent question to determine how best to perform feature selection on these beadarrays. the illumina infinium assay utilizes a pair of probes for each cpg site, one probe for the methylated allele and the other for the unmethylated version. the methylation level is then estimated, based on the measured intensities of this pair of probes. two quantification methods have been proposed to estimate the methylation level: β-values and m-values. while the β-value measures the percentage of methylation at the given cpg site and is the manufacturer's recommended measure to use, the m-value has recently been proposed as an alternative measure. the m-value is defined by the log ratio of the intensities of the methylated to unmethylated probe and has also been used in the context of other methylation array technologies. the m-value is also an analogue of the quantity which has been widely used in expression microarray analysis, although there are two important differences: in the infinium dnam assay, both probes are always measured in the same colour channel, thus dye bias does not need to be adjusted for, and the two probes are measured in the same sample. a recent study based on a titration experiment compared β and m-values and concluded that m-values, owing to their more homoscedastic nature, were a better measure to use. however, this study was limited to one data set and feature selection was only investigated using fold-changes, while statistics and p-values were not considered. the purpose of our study is therefore two-fold: to provide a comprehensive comparison of the performance of β and m-value measures in the context of popular feature selection tools that use actual statistics to rank features, and to compare the performance of some of the most popular feature selection and classification methods on dnam data from different tissues and correlating with different phenotypes, so as to gain objective insights into how best to perform feature selection and classification in this novel context. to address these goals, we make use of independent sets of human dna methylation data, all generated with the 27k illumina infinium platform, encompassing over samples and representing over million data points. materials and methods data sets the human dna methylation data sets are summarised in table all the data were generated using the illumina infinium humanmethylation beadchip assay that enables the direct measurement of methylation at over individual cytosines at cpg loci located primarily in the promoter regions of unique genes. all data sets, except the tcga lung cancer set, followed the same quality control and normalisation strategy. briefly, non background-corrected data was used and only cpgs with intensity detection p-values less than were selected. inter-array normalisation correcting for beadchip effects and variations in bisulfite conversion efficiency were performed using a linear model framework with explicit adjustment for these factors, but only if these factors were significantly correlated with pca/svd components. the independent 27k dnam data sets used in this paper: ukops - uk ovarian cancer population study: whole blood samples from women with ovarian cancer and age-matched healthy controls; endom - normal endometrium and endometrial cancer; cervx - normal cervix and cervical cancer; ovc - ovarian cancer tissue; t1d - type diabetes: whole blood samples from type- diabetics; bc - normal breast tissue and breast cancer; and lc - normal lung tissue and lung cancer. number of cpg probes passing quality control, total number of samples, tissue type, number of normal/cancer samples, age-range and reference to data are given. na: not available. definition of β-value and m-value in the illumina infinium human methylation 27k beadchip assay, bisulphite converted dna is amplified, fragmented and hybridised to the beadchip arrays, with each chip accommodating samples as designated by sentrix positions a-l. each interrogated locus is represented by specific oligomers linked to two bead types, with one representing the sequence for methylated dna and the other for unmethylated dna. the methylation status can be measured either by β-values or m-values. for each specific cpg site, the β-value is calculated from the intensity of the methy and unmethy alleles, as the ratio of the fluorescent signals: β=maxmethy,0/maxmethy,0+maxunmethy,0+ the β-value is a continuous variable between and representing the ratio of the methylated allele to the combined locus intensity. in contrast, the dna methylation m-value is calculated as the log ratio of the intensities of the methylated probe to the unmethylated probe: m=log2maxmethy,0+1/maxunmethy,0+ the m-value is therefore a continuous variable which can in principle take on any value on the real line. since most interrogated cpg sites have intensities larger than, the offset values have only a small effect on both measurements. hence, the relationship between the β and m-values is approximately logistic: β=2∧m/2∧m+1;m=log2β/1-β. while the β-value has a more intuitive biological interpretation it suffers from severe heteroscedasticity. in contrast, the m-value is not directly interpretable in terms of a methylation percentage, but is more statistically valid for analysis of differential methylation levels owing to its more homoscedastic nature. evaluation of feature selection methods we use a multiple training-test set strategy to compare the true positive detection rates of three different feature selection methods: without filtering, filtering based on variability, filtering using supervised principal components. our justification for focusing on vf and spca is that both methods have been very popular and effective in the context of gene expression studies. the performance measure we use to compare the different feature selection methods is the estimated positive predictive value. that is, for each method we select features significantly associated with a phenotype of interest using a training set, and these features are then evaluated in the test-set. those that are also significant in the test set are called true positives, so the ppv measures the probability that a feature declared to be positive in the training set is a true positive. we point out that strictly speaking a true positive can only be a feature which has been shown to be positive using a gold-standard technology. however, in the absence of such gold-standard data, the training/test set strategy advocated here provides a suitable evaluation framework, since a feature that cross-validates is more likely to represent a true positive than one that does not. moreover, this strategy has already been successfully used in previous work, see e.g.. thus, for each of the data sets in this study, independent training and test sets were randomly generated ensuring similar demographics for the phenotype of interest in training and test set. for each training set, a supervised analysis is used to evaluate the association between each cpg methylation profile and the phenotype of interest. in the case of a continuous phenotype the association is carried out under a linear regression model framework, while for a binary phenotype we use t-statistics. test statistic, p-value, q-value, absolute difference in means, and variance of each cpg are recorded. variance-filtering is a feature selection method that filters the cpgs based on their variability before the supervised analysis is performed. the steps are i. in each training set, select the most variable cpgs; ii. rank the most variable cpgs according to the significance of their p-values and retain the top ranked cpgs.. iii. using stepwise evaluation sets, count the number of true positives: i.e. how many cpgs have a statistic p-value < in the corresponding test set while exhibiting the same directional change as that in the training set; iv. for each evaluation set, the positive predictive value is defined as the percentage of cpgs that are true positives. supervised pca is a powerful semi-supervised approach for feature selection and classification, first proposed in. here we adapt it to the feature selection context. the steps are: i. in each training set, perform the supervised analysis and rank all cpgs according to the significance of their p-values; ii. select the top ranked cpgs; iii. run singular value decomposition analysis on the top ranked cpgs and find the principal component most correlated to the phenotype of interest; iv. rank the cpgs according to the absolute value of their corresponding coefficient in the selected principal component and further select the cpgs with the largest absolute weights; v. use the stepwise evaluation sets and count, for each evaluation set size, the number of true positives in the test set using the same criteria as for the vf method; vi. for each evaluation set, compute the ppv as described for vf method. in the without filtering method we perform for each training set the supervised analysis on all ~ cpgs, select the top and then perform the same stepwise evaluation set analysis to estimate ppvs from the p-values in the corresponding test set. the result of applying each of these methods is a matrix of ppv values with rows labelling the evaluation set size and columns labelling the different runs. choice of evaluation set size given that different studies may exhibit widely different operating characteristics, even in relation to the same phenotype of interest, we decided that a more objective comparison of different feature selection methods across studies would be facilitated by choosing evaluation set sizes in each study that result in broadly similar ppv values across studies. hence, evaluation set sizes were chosen so as to minimise differences in ppv between studies, subject to reasonable constraints on the minimum and maximum set sizes. comparison of classification algorithms we use the same multiple training-test partition strategy as used in our feature selection comparison, to evaluate different classification algorithms. thus, within each study, training and test sets were always of the same size and were balanced for the phenotype of interest. we focused on four different powerful classification algorithms, which have been popular in the gene expression field: supervised pca, the lasso algorithm, the elastic net and support vector machines. we used the implementations of these algorithms as provided in the r-packages superpc, glmnet and e in the case of spca we considered classifiers built from up to principal components, and the number of features was optimised in the training sets using internal 10-fold cross-validation. lasso is a special case of the elastic net with the penalty parameter α =. the elastic net itself was run with α =. the additional penalty elnet parameter λ was estimated in the training sets using internal 10-fold cross-validations. the svm was run with a radial kernel with parameter γ = and with error term ε =. these svm parameters were fixed and not optimised, since we wanted to compare the algorithms at the same level of computational complexity. we focused on age as the phenotype of interest since it represents a challenging scenario of small effect sizes and yet it is also well established that age has a significant impact on dna methylation patterns. since age is a continuous variable, svm was run in regression mode, and thus for all methods the predictor is a continuous score. to evaluate concordance of age with the predictor in the test sets we used the c-index. comparison of dimensional reduction algorithms we considered two popular dimensional reduction algorithms: singular value decomposition, and non-negative matrix factorization. in this work we use the svd implementation as computed by the lapack/linpack routine available in r http://www.r-project.org. the application of nmf to dna methylation is justified due to the positivity of beta-valued data. to perform nmf, the "nmf" r package was used. nmf was run using the "brunet" algorithm and initialised using non-negative double svd. hence, using this nndsvd initialisation nmf yields the same solution under repeated runs and allows for a direct comparison to svd. the evaluation of the algorithms was performed by comparing the correlations of the inferred components with the phenotypes of interest, specifically we derived and compared the corresponding r values under a linear regression model. in the case of nmf, we needed to specify the number of components to infer. this number was estimated from the svd analysis by comparing the spectral eigenvalues to the corresponding ones obtained by randomly permuting elements in the data matrix. this number was also the number of svd components used when evaluating correlations with the phenotype of interest.  the signal to noise and signal strength landscape of dna methylation studies one of the aims of this study is to compare different feature selection and quantification methods in relation to important study-specific parameters. one of these parameters is the typical effect size of the study. in the case of a binary phenotype, the effect size of a cpg is defined, loosely speaking, as the ratio of the difference in means between the two phenotypes to a pooled standard deviation. in other words, it can be thought of as a signal to noise ratio, where the noise captures both biological and technical variation. another important parameter to take into consideration in any association analysis is the number of features with an effect size larger than some significance threshold. thus, in this study we ask if the effect size of cpgs associated with a phenotype of interest and their number have an impact on the performance of the different feature selection methods and if this depends on the methylation measure used. in order to consider a wide range of different effect sizes and signal strengths, we considered two main phenotypes: cancer/normal status and age; and two different tissue types: epithelial tissue and whole blood. the corresponding 27k data sets used in this study are summarised in table we verified that these studies exhibited a wide range of different effect sizes and signal strengths, depending largely on tissue type. for instance, we observed, as expected, that cancer diagnostic dnam markers in epithelial tissue have large snrs and signal strengths, in contrast to cancer diagnostic markers in blood, which, while numerous, were characterised by much smaller effect sizes. this fits in with the expectation that cancer associated changes in whole blood reflect mostly changes in the underlying blood cell type composition. in stark contrast to diagnostic markers, we observed that age has a relatively much weaker impact on dnam patterns, and therefore estimated effect sizes were also smaller. the signal strength of age-associated cpgs also varied depending on the tissue in which they were measured. the signal strength and the average effect size of the top ranked cpgs for diagnosis and age-settings and for epithelial and blood tissue. beta-values were used to estimate effect sizes. in the case of age, samples were divided into two groups by the median. m-values significantly outperform β values but only in the small sample size limit a previous study based on a titration experiment advocated the use of m-values over β values, since m-values exhibit less heteroscedasticity. however, this study only used fold-changes to rank features. it is therefore unclear if results differ had the features been ranked using statistics. thus, we asked if using β-values or m-values to quantify methylation has a differential impact on the statistics of differential methylation and on the positive predictive value. first, we verified that in our data sets, β values were highly heteroscedastic, with very low variability at the extremes of the beta-scale, while m-values were approximately homoscedastic. using t-statistics and q-values to estimate the false discovery rate at different p-value thresholds and using all available samples, we observed however that the difference in fdrs between β and m-values within each study was very minor. next, we compared the ppvs obtained using either β or m-values using a multiple 50% training 50% test set strategy where features were ranked according to their t-statistics. because for a given methylation measure we observed substantial differences in fdr rates between studies, in some cases even between studies that looked at the same phenotype of interest and tissue type, we decided to choose evaluation set sizes specific to each study to ensure that ppv values were as similar across studies as possible. in most studies there was no appreciable difference in the ppv between β and m- values and the absolute magnitude of the t-test statistics of the top ranked cpgs were also very similar. to study the dependency on sample size we considered again a multiple training-test set strategy but this time using training and test sets where the number of samples of each phenotype was low. since ordinary t-statistics can not be used in this small sample size limit, we used instead a popular empirical bayes approach to rank features according to a regularized t-statistic. in this small sample size limit, we observed that the ppv was significantly improved when using m-values, in some instances by at least 10%. interestingly, using 50% training/test partitions, corresponding to the largest possible sample sizes, the ppv for m-values was still higher than for β-values, but only by at most 2%. however, for larger sample sizes, ppv values were close to and thus differences between m-values and β-values would naturally be smaller. thus, the improved performance of m-values over β-values is mainly due to the use of regularized t-statistics, since we also observed that in the large sample size limit ordinary t-statistics performed equally well irrespective of m-value or β-value basis. we also observed that in the large sample size limit, ordinary t-statistics performed similarly if not better than regularized t-statistics. thus, given that β-values are more directly interpretable and that in this work we are mostly interested in the performance for relatively large sample sizes, we henceforth focus mainly on β-values. variance filtering can significantly reduce the ppv in the setting of small effect sizes next, we evaluated the effect of filtering cpgs based on variance prior to performing the supervised analysis. variance filtering has been a very popular feature selection tool in gene expression studies, and it has been shown to improve the detection rate of true positives. hence, we asked if the same result holds in the context of dna methylation data. we observed that filtering based on variance improved the ppv in of the cancer/normal studies and in of the studies when age was the phenotype of interest. in the age-setting, vf did worse than wf in two studies by approximately 5% in ppv values, while in the diagnostic setting the worse performance observed in two studies was far less marginal, indicating that in cases where effect sizes and signal strengths are small, filtering based on variance can be counterproductive. thus, we posited that in the case of cancer diagnostic markers in whole blood, that vf would also lead to a decreased ppv performance if larger evaluation sets are used since the larger evaluation sets would include proportionally more false positives. similarly, we would expect that for larger evaluation sets the wf method would allow proportionally more true positives to be captured. confirming this, we observed that wf outperformed vf as the number of top ranked cpgs is increased. this switch in performance between wf and vf as a function of evaluation set size is a consequence of the relatively small effect sizes combined with the relatively large signal strengths of diagnostic markers in whole blood. spca outperforms wf and vf when effect sizes are small supervised principal components has been shown to be a simple yet powerful algorithm for performing feature selection and classification in the gene expression context. we therefore decided to compare spca to the wf and vf methods in terms of feature selection. interestingly, in the diagnostic setting spca did not always outperform vf or wf, although in all data sets it achieved very high ppv values. however, in the more challenging scenario of small effect sizes and signal strengths, as illustrated by using age as the phenotype of interest, we can see that spca achieved substantially better ppv values. thus, principal components can identify true positives more reliably when signal strengths are small. on the other hand, if signal strengths are large, a principal component may not have the plasticity to capture all the important true positives, and thus inevitably will include false positives among the less weighted features within the principal component. we verified that in the diagnostic setting, spca ppv values were substantially higher and comparable to those of wf and vf when using only the top cpgs in the selected pc, demonstrating that at more relaxed feature selection thresholds spca includes numerous false positives signal correlation structure is improved using beta values given that overall feature selection was optimal using spca, we next asked if performance would be similar had we used m-values. remarkably, we observed that spca ppv values were generally higher when evaluated using beta values. this was specially true in the diagnostic setting and was independent of tissue type, indicating that signal strength has a major impact on the inference of biologically relevant principal components. interestingly, therefore, beta-values provide a basis in which the correlative structure of cpg methylation profiles associated with cancer status is better preserved. in contrast, in the m-value basis, it appears that the correlation structure between biologically relevant cpgs is compromised leading to worse modelling of the biological variation. the elastic net and svms outperform spca and lasso next, we turned our attention to the classification task. we considered a total of four different powerful classification algorithms, including spca, lasso regression, the elastic net and support vector machines. we used the same training-test set partition strategy as with our feature selection analysis. age was chosen as the phenotype of interest, since it is now well established that age has an impact on dnam patterns and because age-associated effects are of a relatively small magnitude, thus also providing a more challenging evaluation scenario for the classification algorithms. moreover, age-associated dnam markers may represent cancer risk markers. we considered the same studies as those used for the feature selection analysis in the age setting. across all studies, we observed that either elnet or svm were superior to lasso and spca. the fact that elnet outperformed lasso in all data sets extends previous results obtained on gene expression to the dna methylation context. moreover, we observed that spca always performed optimally in the context of a single component, i.e. adding more principal components to the classifier did not improve performance. interestingly, spca- performed similarly to elnet and svm in the two smaller studies, suggesting that the rigidity imposed by principal components can be of an advantage in this setting. unsupervised modelling of diagnostic and age effects finally, we compared unsupervised algorithms in their ability to model diagnostic and age effects in dna methylation studies. we considered two of the most popular unsupervised dimensional reduction algorithms from the gene expression field: singular value decomposition, and non-negative matrix factorization. given that dna methylation beta-valued data are defined on the compact support and are thus positively valued, it is justified to explore the application of nmf in this context. to objectively compare nmf to svd, nmf was initialised using non-negative double svd, thus allowing us to determine if the positivity constraints of the nmf framework add biological value. for both svd and nmf we computed the r correlations between the inferred components and the phenotype of interest. as a modelling criterion we compared the average r of the two best components, obtained from each algorithm. focusing first on the diagnostic setting, we can clearly observe that nmf outperformed svd/pca in all studies, including the study comparing blood samples from cancers and controls. however, nmf was not conclusively better in the age-setting of much smaller effect sizes. discussion a previous study compared β to m-values in the context of a titration experiment using fold-change values to rank cpgs. the study concluded that m-values, owing to their reduced heteroscedasticity, provided a better basis in which to detect true positives. this result is not unexpected because fold changes only consider the magnitude of differences in mean values, regardless of the underlying biological and technical variability. therefore, using β-values and fold changes, true positive cpgs with low or high mean beta methylation values would not be highly ranked and would result in a significant proportion of false negatives. as we have shown here, using many different data sets and using a multiple training test set strategy to estimate the ppv, t-statistics derived from m-values and β-values are comparable and lead to similar ppv values, but only in the scenario of relatively large sample sizes. in contrast, when sample sizes are small and a regularized t-statistic must be used, m-values provided a significantly better basis for inference, sometimes by as much as 10% improvements in the ppv. that regularized t-statistics perform relatively poorly in the β-value basis is not surprising, since the regularisation involves an empirical bayes approach in which a posterior variance is estimated from estimates specifying the prior, an approach which is known to be more sensitive to the precise distributional properties of the data and is therefore also more sensitive to deviations from normality. therefore, our comprehensive study using statistics to rank features, indicates that the severe heteroscedasticity of β-values does not compromise the reliability of identifying true positives in the case of relatively large sample sizes, but that m-values offer significantly improved inference in the limit of small sample sizes. a second important finding is that variance filtering does not always increase the ppv. indeed, in scenarios where effect sizes are small but where there are numerous associated features, we have seen how variance filtering can compromise the detection of true positives when more relaxed significance thresholds are used. this was illustrated in the context of cancer diagnostic methylation markers in whole blood, where the changes reflect underlying changes in blood cell type composition as demonstrated by us previously. in this setting, effect sizes are small but may involve a potentially large number of genes, indeed potentially all those genes whose expression differentiates blood cell types from each other and whose expression is under epigenetic regulation. since very often there is a desire to perform gene set enrichment analysis on a set of top ranked features, thresholds must be chosen to ensure a relatively small false negative rate. hence, in the scenario described here, variance filtering may be specially counterproductive since it leads to a larger fnr and hence could compromise the detection of enriched go-terms or pathways. it follows that any decision to filter features based on variance must take into consideration the clinical and biological context of the study. our third key result is the improved performance of spca in detecting true positives, specially in the more challenging scenario when effect sizes and signal strengths are small. while spca has been shown to outperform many other feature selection methods in the context of gene expression data, an analogous result was still lacking in the case of dna methylation. our data therefore show that spca is a powerful feature selection algorithm, independently of the molecular profiles considered or their underlying statistical distributions. the improved performance over vf also indicates that the correlative dnam patterns can be exploited to identify true positives more reliably. interestingly, we observed that spca worked better in the context of β-values. this may be surprising as one might be inclined to believe that components of maximal variation are more robustly identified in a basis that promotes variation, i.e. using m-values. however, it is also the case that m-values would aggravate the effect of outliers and therefore distort the true components of maximal biological variation. therefore, illumina's recommendation to use β-values over m-values appears to be further justified by the fact that β-values, by virtue of being bounded, may provide a natural regularization and thus suppress the undesirable effects of potential technical outliers. we also compared various powerful classification algorithms, including spca, elnet, lasso and svms. our results clearly indicate that elnet and svm are superior classification methods in the case where effect sizes are small. given that elnet also provides an automatic means of feature selection, it would appear to be a preferable choice over svm. nevertheless, svm outperformed elnet in three studies, including the two smaller ones. interestingly, we also observed that classification using spca performed optimally in the context of a single component of variation. since multi-component spca uses a training set to estimate an optimal linear combination of principal components that then makes up the predictor, it would appear that the estimated weights specifying the linear combination is not mirrored in the test set, hence why spca- is optimal. it will be interesting to investigate if this result generalizes to other phenotypes where effect sizes may be small. finally, we also compared two of the most popular unsupervised algorithms in their ability to model differences between cancer and normal tissues and differences associated with age. similar to the results obtained on gene expression, we observed that nmf also outperforms svd/pca in the context of dna methylation data. indeed, the meta-genes inferred from the nmf were generally more highly correlated with cancer/normal status across all data sets considered. hence, despite the increased computational complexity of nmf, this algorithm should be used whenever expected effect sizes are large. on the other hand, if effect sizes are small we did not find that the increased computational complexity of nmf offered any advantage over pca. in the context of evaluating feature selection and classification methods, it is important to discuss the impact of potential batch effects. indeed, as shown in many previous papers, known and unknown batch effects can affect a substantial proportion of features in an experiment and lead to biased estimates of statistical significance. in this work we used inter-array normalised data where beadchip effects and variations in bisulfite conversion efficiency were adjusted for using a linear model framework, as these factors were observed to account for significant components of variation in pca analysis. we observed that not correcting for these factors may alter the absolute performance of any given method, yet the relative performance of the different methods was largely robust to whether inter-array normalisation was performed or not. since the inter-array normalisation used only adjusts for known effects, we have not considered here the potential effects of unknown confounders. to address this would require adapting surrogate variable analysis approaches in the training/test set evaluation framework considered here. this has only very recently made possible through an extension of sva, called frozen sva. thus, it will be interesting to further compare these feature selection methods in the context of the present dna methylation data sets, and thus further assess the potential impact of batch effects on feature selection and classification. while the insights obtained in this study have been derived from one particular array, it is likely that similar insights would apply to the recently released scaled up 450k infinium beadchip. we should point out however, that one of the key differences between the 27k and 450k platforms is that the probes on the 450k beadarray come in two different designs and are thus characterised by widely different statistical distributions. thus, the insights obtained here are likely to apply only to the data restricted to probes of one particular design. in any case, the much higher density of the 450k array, means that feature selection in 450k data will be even more critical.  the optimal choice of methylation measure and feature selection method in the context of 27k dna methylation data is dependent on the sample size, expected effect sizes and signal strength of the specific study. while spca is a powerful feature selection tool in dna methylation studies where effect sizes and signal strengths are small, the elastic net and svms generally outperform spca in the context of classification. these insights are important for anyone embarking on large scale dna methylation profiling using 27k or 450k beadarrays. abbreviations dnam: dna methylation; pca: principal components analysis; svd: singular value decomposition; spca: supervised principal components; elnet: the elastic net; svm: support vector machines; nmf: non-negative matrix factorisation. competing interests the authors declare that they have no competing interests. authors' contributions jz and aet performed the statistical analyses. mw contributed samples and obtained funding. aet wrote the manuscript with contributions from jz. all authors read and approved the final manuscript supplementary material additional file 1 typical scatterplots of mean methylation vs standard deviation in methylation using beta and m values. green line denotes a mean loess smoother. a) the blood samples from the healthy controls in the ukops study. b) the blood samples from the t1d study. click here for file additional file 2 false discovery rate estimation table at different thresholds of top cancer diagnostic associated cpgs obtained from the supervised analysis in the ukops, endom, and cervx data sets. click here for file additional file 3 the evaluation set size is plotted against the mean ppv for the different studies. the mean ppv represents an average over the training-test partitions and over the feature selection methods. a) diagnostic setting, b) age setting. the selected evaluation set sizes in each study and for each phenotype of interest are marked in red. note that in the case of age, differences in mean ppv between studies could not be minimized because of the additional constraint of a reasonable minimum set size. click here for file additional file 4 a-b) comparison of positive predictive values obtained from beta and m-values using the without-filtering method and using training-test set partitions. a) diagnosis setting in ukops, endom, cervx, bc and lc with evaluation set sizes of top and diagnosis-associated-cpgs, respectively; and b) age setting in t1d, ovc, cervx, ukops, endom, and bc with the top and age-associated-cpgs as evaluation sets, respectively. c-d) comparison of absolute t-statistics obtained from beta and m-values. c) diagnostic setting in ukops, endom, cervx, bc and lc; and d) age setting in t1d, ovc, cervx, ukops, endom, and bc. all p-values shown are from a wilcoxon-rank sum test. click here for file additional file 5 mean ppv values between spca, wf and vf in the diagnostic setting but using only the top ranked cpgs as evaluation sets. click here for file acknowledgements this work was funded by the eve appeal and a grant from the uclh/ucl comprehensive biomedical research centre project and undertaken at uclh/ucl, which received a proportion of its funding from the department of health nihr biomedical research centres funding scheme. this work was also supported by the european network translational research in gynaecological oncology of the european society of gynaecological oncology. we thank allison jones, usha menon, helga salvesen and heidi fiegl for providing data. aet is supported by a heller research fellowship.
our results highlight the importance of tailoring the feature selection and classification methodology to the sample size and biological context of the dna methylation study. we show that for relatively large sample sizes feature selection using test statistics is similar for m and β-values, but that in the limit of small sample sizes, m-values allow more reliable identification of true positives. finally, in unsupervised modelling of cancer diagnosis, we find that non-negative matrix factorisation clearly outperforms principal components analysis. we also demonstrate that classification using the elastic net and support vector machine clearly outperforms competing methods like lasso and spca. the elastic net emerges as a powerful classification algorithm for large-scale dna methylation studies, while nmf does well in the unsupervised context. using a total of large illumina infinium 27k methylation data sets, encompassing over samples from a wide range of tissues, we here provide an evaluation of popular feature selection, dimensional reduction and classification methods on dna methylation data. specifically, we find that variance filtering improves the detection of true positives in studies with large effect sizes, but that it may lead to worse performance in studies with smaller yet significant effect sizes. while feature selection and classification methods have been comprehensively explored in the context of gene expression data, relatively little is known as to how best to perform feature selection or classification in the context of illumina infinium methylation data. we also show that the effect of variance filtering on feature selection is study-specific and dependent on the phenotype of interest and tissue type profiled. in contrast, supervised principal components improves the statistical power, especially in studies with small effect sizes. given the rising importance of epigenomics in cancer and other complex genetic diseases, and in view of the upcoming epigenome wide association studies, it is critical to identify the statistical methods that offer improved inference in this novel context. specifically, we evaluate the effects of variance filtering, supervised principal components and the choice of dna methylation quantification measure on downstream statistical inference.
 the 27k illumina infinium methylation beadchip is a popular high-throughput technology that allows the methylation state of over cpgs to be assayed. while feature selection and classification methods have been comprehensively explored in the context of gene expression data, relatively little is known as to how best to perform feature selection or classification in the context of illumina infinium methylation data. given the rising importance of epigenomics in cancer and other complex genetic diseases, and in view of the upcoming epigenome wide association studies, it is critical to identify the statistical methods that offer improved inference in this novel context.  using a total of large illumina infinium 27k methylation data sets, encompassing over samples from a wide range of tissues, we here provide an evaluation of popular feature selection, dimensional reduction and classification methods on dna methylation data. specifically, we evaluate the effects of variance filtering, supervised principal components and the choice of dna methylation quantification measure on downstream statistical inference. we show that for relatively large sample sizes feature selection using test statistics is similar for m and β-values, but that in the limit of small sample sizes, m-values allow more reliable identification of true positives. we also show that the effect of variance filtering on feature selection is study-specific and dependent on the phenotype of interest and tissue type profiled. specifically, we find that variance filtering improves the detection of true positives in studies with large effect sizes, but that it may lead to worse performance in studies with smaller yet significant effect sizes. in contrast, supervised principal components improves the statistical power, especially in studies with small effect sizes. we also demonstrate that classification using the elastic net and support vector machine clearly outperforms competing methods like lasso and spca. finally, in unsupervised modelling of cancer diagnosis, we find that non-negative matrix factorisation clearly outperforms principal components analysis.  our results highlight the importance of tailoring the feature selection and classification methodology to the sample size and biological context of the dna methylation study. the elastic net emerges as a powerful classification algorithm for large-scale dna methylation studies, while nmf does well in the unsupervised context. the insights presented here will be useful to any study embarking on large-scale dna methylation profiling using illumina infinium beadarrays. dna methylationclassificationfeature selectionbeadarrays 