# An Efficient Approach to Screening Epigenome-Wide Data

Authors: Meredith A. Ray; Xin Tong; Gabrielle A. Lockett; Hongmei Zhang; Wilfried J. J. Karmaus
Journal: BioMed Research International
 (2016)

---

## Abstract

Screening cytosine-phosphate-guanine dinucleotide (CpG) DNA methylation sites in association with some covariate(s) is desired due to high dimensionality. We incorporate surrogate variable analyses (SVAs) into (ordinary or robust) linear regressions and utilize training and testing samples for nested validation to screen CpG sites. SVA is to account for variations in the methylation not explained by the specified covariate(s) and adjust for confounding effects. To make it easier to users, this screening method is built into a user-friendly R package,ttScreening, with efficient algorithms implemented. Various simulations were implemented to examine the robustness and sensitivity of the method compared to the classical approaches controlling for multiple testing: the false discovery rates-based (FDR-based) and the Bonferroni-based methods. The proposed approach in general performs better and has the potential to control both types I and II errors. We appliedttScreening to 383,998 CpG sites in association with maternal smoking, one of the leading factors for cancer risk.


---

## Body

## 1. Background

Due to its high throughput, accuracy, small sample requirement, and acceptable cost, the Illumina Infinium HumanMethylation450 BeadChip has been widely used to analyze deoxyribonucleic acid (DNA) methylation profiles in epigenetic studies that target various types of cancer. In particular, the illumina infinium assay utilizes a pair of probes (a methylated probe and an unmethylated probe) to measure the intensities of methylated and unmethylated alleles at the interrogated cytosine-phosphate-guanine dinucleotide (CpG) sites [1]. Two measures of DNA methylation are usually used: beta-values and M-values. A beta-value is the ratio of signal from the methylated probe relative to the sum of both methylated and unmethylated probes. Beta-values are in the range of (
0,1
) with 0 being completely unmethylated and 1 being fully methylated. M-values are log
⁡
2 ratio of intensities for methylated and unmethylated probes and range from (
-
∞
,
+
∞
) [2, 3]. M-values are used more often in appreciation of its wide data range and variance homogeneity compared to beta-values.Given the feature of high dimensionality of high-throughput methylation data, when performing designed and possibly complicated statistical analyses, it is wise to target potentially important CpG sites, for instance, CpG sites potentially associated with single-nucleotide polymorphisms (SNPs) and/or other covariate(s) of interest. Otherwise, the statistical power will be substantially lost. There is evidence that methylation is affected by genetic and some known factors such as smoking [4, 5]. Screening CpG sites have become overwhelmingly important across multiple health fields of study such as cancer research, genetic diseases, and epigenetic research.Common methods for screening CpG sites assume some relationship of theM-values in association with SNPs or some other genetic or environmental factors conditional on the assumption of linearity with some post hoc adjustment for multiple comparisons. The advantage to this method is the flexibility of incorporating additional covariates and their interactions. The primary limitation lies in controlling for multiple testing. Two popular adjustment methods are the Bonferroni-based method [6, 7] and the Benjamini-Hochberg method for controlling the false discovery rate (FDR) [8, 9]. These methods alter the p-value or critical value to control for type I error. Bonferroni correction is the most conservative by dividing the linear regression p-value, respective of the regression term of interest, by the total number of comparisons (
m
) or CpG sites in this case such that those adjusted p-values above the significance level are rejected. The FDR method first orders the p-values, P
(
k
) for k
∈
1
⋯
m, such that lower ordered p-values that are less than or equal to k
/
m
×
α are rejected [9]. It follows that the conservative Bonferroni-based method cannot control for type II error while the FDR-based method cannot control for type I error.There are other potential issues that arise when dealing with DNA methylation. It is possible that the variation in methylation cannot be fully explained by the known covariates and there exist latent factors that confound with these known covariates [5]. To improve the screening quality, it is thus important to account for variations introduced by other unknown factors. Furthermore, CpG sites screened from one data set may not be consistent with those from another data set which directly affects the type I error rate and leads to a loss of power. It is thus equally important to improve the reproducibility of the selected CpG sites.In this paper, we propose a novel collaboration of existing statistical techniques to screen genome-wide methylation. It takes unknown factor effects into account and achieves better reproducibility. The method has the ability to control for both types I and II error while adjusting for covariates as well as latent variables. The proposed screening method incorporates surrogate variable analysis [5], which identifies unknown latent variables, in conjunction with a training and testing approach [10] across CpG methylation sites linearly associated with covariates of interest, including the identified surrogate variables. Independently, each method is well established for different purposes. We mingle these methods and form, compared to existing methods, an improved and more efficient process to screen (filter) informative DNA methylation sites. In addition, this proposed method has been built into an efficient and user-friendly R package. In the following sections, further description and details of the proposed method can be found in Section 2, simulation studies and a real data application are included in Section 3, and we summarize the approach in Section 4.


## 2. Materials and Methods

The proposed screening procedure is built for analyzing the associations between methylation data (M-values), or some high-dimensional data, and covariates of interest and their potential interactions. It consists of two consecutive components and surrogate variable analysis followed by a series of regressions while controlling for multiple testing. Surrogate variable analysis (SVA) aims to identify and estimate latent factors or surrogate variables (SVs) that potentially affect the association between known factors and the response variable, for example, SNPs (known factor) and DNA methylation (response variable) [5]. Including the estimated surrogate variables into the screening process has the potential to reduce unexplained variations, adjust for confounding effects, and consequently improve the accuracy of screening in terms of important variable identification [5]. In the context of DNA methylation, inclusion of the surrogate variables explains the variation in DNA methylation not explained by the covariates currently under consideration. This implies that the identified surrogate variables can be further used to identify important factors (markers) showing large contribution to the variation in the response variable explained by the surrogate variables [11]. These surrogate variables along with other variable(s) of interest can be included in regression analysis as independent variables. After SVA, we then begin the screening process with regressions and adjust for multiple comparisons.As mentioned in the introduction, several methods exist to adjust for multiple testing but lack the ability to control for both types I and II errors. We have elected to implement a method that will control for both while simultaneously helping the issue of reproducibility. We let randomly chosen training and testing samples estimate and test the effects of the primary covariate(s), termed the TT method. General ideas of this approach are discussed in Dobbin and Simon and Faraggi and Simon [10, 12]. This method follows the concept of cross-validation. It has been shown that the implementation of the training and testing technique can provide a better control of type I error rate [10, 12]. In the next two subsections, we provide detailed steps and options available to both the surrogate variable analysis component and the TT method.


### 2.1. Identifying Surrogate Variables

Surrogate variables are inferred prior to screening using an algorithmic method developed by Leek and Storey (2007) called surrogate variable analysis [5]. Following the descriptions in Leek and Storey (2007), these SVs are developed by removing the amount of methylation or signal due to the variable(s) of interest and then decomposing the remaining residuals to identify an orthogonal basis of singular vectors that can be reproduced. These vectors are further examined for significant variation to form surrogate variables. Leek and Storey built an R package to perform the surrogate variable analyses (SVAs). The first step in SVA is to identify the number of surrogate variables based on the data using one of two methods, “be” or “leek” as noted in the package. The “be” method, being the default choice according Leek and Storey’s R package, is based on a permutation procedure originally proposed by Buja and Eyuboglu in 1992 [13], while the “leek” method provides an interface to the asymptotic approach proposed by Leek in 2011 [14]. Once the number of surrogate variables is calculated, they are then estimated using one of three algorithms, the iteratively reweighted (“irw”), supervised (“supervised”), and the two-step (“two-step”) method. Iteratively reweighted method is for empirical estimation of control probes, supervised method is for when control probes are known, and the two-step method is for general estimation of surrogate variables [15]. We elected to implement the “two-step” method following Leek and Storey [5]. Conditional on the data, a number of latent unknown variables have been identified, estimated, and will now be incorporated into the regression in association with DNA methylation as additional covariates.


### 2.2. The Training and Testing (TT) Screening Method

After surrogate variable analysis is complete, the TT screening method then begins as an iterative process of randomly sampling the training and testing data by a specified proportion. By default, 2/3 of the data will be included in the training data set, which is suggested in Dobbin and Simon [10] to maximize statistical power. Linear regressions are first applied to the training data to calculate the p values for the association between the CpG site and the covariates, including SVs, using either ordinary least squares (this is a default choice in our R package) or robust regression. Robust regression is a type of linear regression that allows for more relaxed assumptions about normality and presence of outliers in the data. A CpG site is included as a candidate if the covariate(s) of interest is statistically significant according to a prespecified significance level, for example, 0.05. In our designed R package, we give user the flexibility to define which term (covariate) is used to decide the selection of CpG sites. For example, suppose the defined right-hand side of the regression is x
1
+
x
2
+
x
1
×
x
2, where x
1
×
x
2 denotes the interaction of x
1 and x
2. If the decision of selecting a CpG site is based on one single term, for example, the significance of the interaction effect, the p value for the interaction term would be used to test statistical significance against the prespecified significance level of 0.05.The process continues with these candidate CpG sites being further tested using the remaining subjects (testing data set) with linear regressions. For one pair of training and testing data sets, a candidate CpG site is deemed as being important if the significance still holds in the testing data. The significance level for the testing data by default is set at the same level as for the training data, 0.05. This screening process will be repeatedi times (i.e., iterations); at each iteration, a training and testing data set will be randomly selected. After one iteration, a pool of candidate CpG sites will be selected. This process is continued for total i iterations. We summarize this screening process in Figure 1. Across all i iterations, CpG sites selected in at least m iterations will be included in the final pool of potentially important CpG sites; that is, the cutoff percentage of selection is m
/
i
×
100
%. Final estimates of the associations and statistical significance for the selected CpG sites are inferred by use of the complete data via the same analytical methods (i.e., linear regression including surrogate variables previously estimated from complete data) as in the training and testing process. A cutoff percentage of 50 (m
=
50 across 100 total iterations (i
=
100)) was used to determine the final pool of potentially important CpG sites. Suggestions on the determination of this predefined value, m, are discussed later in Section 3.1.Figure 1
A conceptual representation of the training and testing algorithm.The user-friendly R package,ttScreening, is available from the Comprehensive R Archive Network (CRAN) at https://cran.r-project.org/web/packages/ttScreening/index.html, which implements the proposed screening procedure discussed above. ThisttScreening package also provides access to other screening methods: FDR and Bonferroni methods. Various options, such as type of linear regression and surrogate variable estimation method, are available for the user to specify while other options are data specific and will need to be defined by the user. However, the package does provide acceptable defaults values for those options that are not data specific. A list of package options along with descriptions are available in the package manual at https://cran.r-project.org/web/packages/ttScreening/ttScreening.pdf.


## 2.1. Identifying Surrogate Variables

Surrogate variables are inferred prior to screening using an algorithmic method developed by Leek and Storey (2007) called surrogate variable analysis [5]. Following the descriptions in Leek and Storey (2007), these SVs are developed by removing the amount of methylation or signal due to the variable(s) of interest and then decomposing the remaining residuals to identify an orthogonal basis of singular vectors that can be reproduced. These vectors are further examined for significant variation to form surrogate variables. Leek and Storey built an R package to perform the surrogate variable analyses (SVAs). The first step in SVA is to identify the number of surrogate variables based on the data using one of two methods, “be” or “leek” as noted in the package. The “be” method, being the default choice according Leek and Storey’s R package, is based on a permutation procedure originally proposed by Buja and Eyuboglu in 1992 [13], while the “leek” method provides an interface to the asymptotic approach proposed by Leek in 2011 [14]. Once the number of surrogate variables is calculated, they are then estimated using one of three algorithms, the iteratively reweighted (“irw”), supervised (“supervised”), and the two-step (“two-step”) method. Iteratively reweighted method is for empirical estimation of control probes, supervised method is for when control probes are known, and the two-step method is for general estimation of surrogate variables [15]. We elected to implement the “two-step” method following Leek and Storey [5]. Conditional on the data, a number of latent unknown variables have been identified, estimated, and will now be incorporated into the regression in association with DNA methylation as additional covariates.


## 2.2. The Training and Testing (TT) Screening Method

After surrogate variable analysis is complete, the TT screening method then begins as an iterative process of randomly sampling the training and testing data by a specified proportion. By default, 2/3 of the data will be included in the training data set, which is suggested in Dobbin and Simon [10] to maximize statistical power. Linear regressions are first applied to the training data to calculate the p values for the association between the CpG site and the covariates, including SVs, using either ordinary least squares (this is a default choice in our R package) or robust regression. Robust regression is a type of linear regression that allows for more relaxed assumptions about normality and presence of outliers in the data. A CpG site is included as a candidate if the covariate(s) of interest is statistically significant according to a prespecified significance level, for example, 0.05. In our designed R package, we give user the flexibility to define which term (covariate) is used to decide the selection of CpG sites. For example, suppose the defined right-hand side of the regression is x
1
+
x
2
+
x
1
×
x
2, where x
1
×
x
2 denotes the interaction of x
1 and x
2. If the decision of selecting a CpG site is based on one single term, for example, the significance of the interaction effect, the p value for the interaction term would be used to test statistical significance against the prespecified significance level of 0.05.The process continues with these candidate CpG sites being further tested using the remaining subjects (testing data set) with linear regressions. For one pair of training and testing data sets, a candidate CpG site is deemed as being important if the significance still holds in the testing data. The significance level for the testing data by default is set at the same level as for the training data, 0.05. This screening process will be repeatedi times (i.e., iterations); at each iteration, a training and testing data set will be randomly selected. After one iteration, a pool of candidate CpG sites will be selected. This process is continued for total i iterations. We summarize this screening process in Figure 1. Across all i iterations, CpG sites selected in at least m iterations will be included in the final pool of potentially important CpG sites; that is, the cutoff percentage of selection is m
/
i
×
100
%. Final estimates of the associations and statistical significance for the selected CpG sites are inferred by use of the complete data via the same analytical methods (i.e., linear regression including surrogate variables previously estimated from complete data) as in the training and testing process. A cutoff percentage of 50 (m
=
50 across 100 total iterations (i
=
100)) was used to determine the final pool of potentially important CpG sites. Suggestions on the determination of this predefined value, m, are discussed later in Section 3.1.Figure 1
A conceptual representation of the training and testing algorithm.The user-friendly R package,ttScreening, is available from the Comprehensive R Archive Network (CRAN) at https://cran.r-project.org/web/packages/ttScreening/index.html, which implements the proposed screening procedure discussed above. ThisttScreening package also provides access to other screening methods: FDR and Bonferroni methods. Various options, such as type of linear regression and surrogate variable estimation method, are available for the user to specify while other options are data specific and will need to be defined by the user. However, the package does provide acceptable defaults values for those options that are not data specific. A list of package options along with descriptions are available in the package manual at https://cran.r-project.org/web/packages/ttScreening/ttScreening.pdf.


## 3. Results and Discussion

Simulations are used to demonstrate and assess the TT screening method in comparison with the FDR- and Bonferroni-based methods using thettScreening R package. These are followed up by an application to a data set of 383,998 CpG sites and their association with maternal smoking.


### 3.1. Simulations

Simulation Scenarios. In total, 2,000 CpG sites across n
=
600 subjects were simulated. In these 2,000 CpG sites, k sites were assumed to be important. Different settings of k were considered, k
=
10, 100, 200, and 400. Among the k important sites, DNA methylation at 90
% of the k CpG sites was associated with two variables x
1 and x
2, their interaction, and 5 unobservable independent uncorrelated variables, and the remaining 10
% of the k sites were associated with x
1, the interaction of x
1 and x
2, and the 5 unobservable independent uncorrelated variables. Variable x
1 is generated from normal distribution with mean 1 and variance 1, and x
2 is a four-level categorical variable generated from multinomial with parameters n and π
=
{
0.15,0.25,0.25,0.35
}. The remaining CpG sites (2000
-
k) were only associated with the 5 unobservable variables and were deemed as unimportant ones. Linear regressions were applied to simulate the data. To assess the robustness of each method (TT screening and FDR- and Bonferroni-based methods), we considered two types of random error in the regressions, one following normal distribution with mean 0 and variance σ
2
=
1.5 and the other χ
2 distribution with a degree freedom of 1. Combining the choices on k and the distributions of random error, in total, we have 8 settings. For each setting, we generated 100 Monte Carlo (MC) replicates. The results presented below include means of the number of incorrect selections (rounded to the nearest integer), estimates of sensitivity, and estimates of specificity across the 100 MC replicates. The number of incorrect selections refers to the number of CpG sites misidentified (refers to both false positive and false negative CpG sites) out of the 2,000 CpG sites.


### 3.2. Results

Variablesx
1 or x
2 would be selected if their interaction effect was statistically significant. All package options were chosen as the default values with training and testing significance levels both set to 0.05. The screening results (Table 1) indicate that, in general, the sensitivity from the FDR-based method is comparable to that from the TT screening method but its specificity is lower when the number of important variables is not sparse; for example, k
=
200. Compared to the Bonferroni-based method, the TT method in general gave better sensitivity and comparable specificity. These were as expected, as the Bonferroni-based method lacks the ability to control type II error while the FDR-based method cannot control well type I error [16]. The TT screening method, on the other hand, has the potential to control both types I and II errors. This is reflected by the results that the TT screening method overall produced the smallest number of incorrectly identified variables, a statistic incorporating information from both sensitivity and specificity. We also performed the screening using robust regressions and similar results were obtained (not shown). These findings are invariant to the distribution pattern of random errors, normally, or skew distributed.Table 1
Simulation results for selectingk important variables among 2,000 candidates with surrogate variables included.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
100
# incorrect
0
1
1
3
8
2
Sensitivity
0.981
0.991
0.994
0.968
0.999
0.991
Specificity
1
1
1
1
0.996
1
k
=
200
k
=
400
# incorrect
13
32
7
78
213
46
Sensitivity
0.936
0.998
0.982
0.814
0.996
0.931
Specificity
1
0.982
0.998
0.998
0.868
0.989
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
100
# incorrect
1
1
1
16
9
7
Sensitivity
0.892
0.958
0.963
0.837
0.982
0.941
Specificity
1
1
1
1
0.996
1
k
=
200
k
=
400
# incorrect
46
29
21
169
158
95
Sensitivity
0.771
0.981
0.906
0.582
0.969
0.790
Specificity
1
0.986
0.999
0.999
0.909
0.993
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.Recall that the default value of the cutoff percentage required for a CpG site to be treated as an informative site across all iterations was 50%, and the significance levels for both the training and testing data set were at 0.05. To further evaluate how the choices of this cutoff percentage and the significance levels in the training and testing steps influence the screening results, we chose a set of the cutoff percentage values ranging from 30 to 90 and set the training significance level to 0.1 instead of default value 0.05. All other settings were kept the same. As seen in Figure2, across the different number of important CpG sites and sample sizes, overall taking the cutoff percentage close to 50% works the best with significance level of 0.05 for both the training and testing steps. If a higher significance level is chosen in the testing step, then a higher value for cutoff percentage should be used.Figure 2
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites. For the TT screening method, two sets of significance levels are considered: (1) 0.05 for training data and 0.1 for testing data; (2) 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000The above examination on cutoff percentage is with sample sizen
=
600. To assess whether and how sample size influences the choice of cutoff percentage, we repeated the above analyses for different sample sizes (and for each sample size, 100 MC replicates were generated), n
=
200, 400, and 800 with all options inttScreening() set at default values except for cutoff percentage. The results (Figures 3–5) indicated that if the number of important variables is expected to be sparse (e.g., 0.5% of the total number variables), the default cutoff value (50%) is a reasonable choice regardless of the sample size. On the other hand, if important variables are not sparse, for example, at least 5% of the total number of variables, then the cutoff value tends to be influenced by the sample size only if sample size is not large, for example, ≤400. In this case, the smaller the sample size is, the lower the cutoff percentage should be taken. This is as expected; when the sample size is smaller, the probability of true positives will be lower. We thus do not expect a higher proportion of correct selection and consequently a lower cutoff percentage is desired. From our simulations, we recommend 20% cutoff if both the following two conditions are not satisfied: (1) important variables are not expected to be sparse (which rarely happens in high-throughput and high-dimensional data) and (2) we have small samples compared to the number of candidate CpG sites. However, if sample size is large, the default cutoff 50% is still recommended. The selection results for sample sizes n
=
200, 400, and 800 following these recommendations are included in the appendix (Tables 3, 4, and 5).Figure 3
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites across 200 subjects. For the TT screening method, significance levels considered are 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000Figure 4
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites across 400 subjects. For the TT screening method, significance levels considered are 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000Figure 5
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites across 800 subjects. For the TT screening method, significance levels considered are 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000To demonstrate the benefit of using the identified surrogate variables in screening, we reanalyzed the simulated data (forn
=
600) but with the surrogate variable analysis excluded (this option is available in the package). The relative patterns of statistics between different methods were intact compared to the results when surrogate variables were included in the screening process. However, when surrogate variables were not considered, a drastic drop in the mean sensitivity was observed across all the settings (Table 2). The sensitivity measures across all the 100 replicates ranged from 0 to 30%, implying that all methods had trouble deciphering CpG sites that are truly important. However, instead of completely excluding surrogate variables, what if we include a set of surrogate variables with each surrogate variable explaining quite a small portion of variation (i.e., less informative)? To examine this, we extracted surrogate variables that are less informative and reanalyzed the data. Similar findings as shown in Table 2 were concluded (Table 6), indicating the importance of including informative surrogate variables in order to substantially improve the quality of selection.Table 2
Simulation results for selectingk important variables among 2,000 candidates with surrogate variables excluded.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
100
# incorrect
10
51
9
99
114
90
Sensitivity
0.011
0.054
0.100
0.014
0.167
0.102
Specificity
1
0.979
1
1
0.984
1
k
=
200
k
=
400
# incorrect
197
182
180
394
310
359
Sensitivity
0.015
0.232
0.102
0.015
0.288
0.104
Specificity
1
0.984
1
1
0.984
1
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
100
# incorrect
10
38
9
98
103
91
Sensitivity
0.018
0.053
0.095
0.015
0.158
0.096
Specificity
1
0.986
1
1
0.99
1
k
=
200
k
=
400
# incorrect
197
176
181
394
314
363
Sensitivity
0.014
0.217
0.094
0.014
0.273
0.095
Specificity
1
0.989
1
1
0.986
1
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.Table 3
Simulation results for selectingk important variables among 2,000 candidates including surrogate variables across 200 subjects.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
10
0
†
# incorrect
7
6
6
77
43
41
Sensitivity
0.292
0.423
0.475
0.227
0.605
0.641
Specificity
1
1
1
1
0.998
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
165
87
93
361
214
244
Sensitivity
0.175
0.616
0.571
0.099
0.565
0.432
Specificity
1
0.994
0.996
1
0.975
0.989
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
10
0
†
# incorrect
9
8
7
88
62
55
Sensitivity
0.126
0.193
0.307
0.124
0.408
0.498
Specificity
1
1
1
1
0.998
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
182
123
119
380
274
288
Sensitivity
0.093
0.422
0.435
0.051
0.377
0.314
Specificity
1
0.996
0.997
1
0.985
0.992
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.†Cutoff percentage of 20%.Table 4
Simulation results for selectingk important variables among 2,000 candidates including surrogate variables across 400 subjects.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
10
0
†
# incorrect
2
1
1
24
10
9
Sensitivity
0.851
0.949
0.941
0.761
0.968
0.965
Specificity
1
1
1
1
0.997
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
64
29
21
210
141
93
Sensitivity
0.682
0.968
0.947
0.479
0.948
0.864
Specificity
1
0.988
0.994
0.999
0.925
0.976
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
10
0
†
# incorrect
4
3
3
47
19
18
Sensitivity
0.612
0.767
0.786
0.531
0.878
0.883
Specificity
1
1
1
1
0.997
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
111
42
42
289
146
141
Sensitivity
0.445
0.88
0.837
0.281
0.841
0.717
Specificity
1
0.99
0.995
0.999
0.949
0.983
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.†Cutoff percentage of 20%.Table 5
Simulation results for selectingk important variables among 2,000 candidates including surrogate variables across 800 subjects.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
100
# incorrect
2
44
5
26
309
36
Sensitivity
0.991
1
0.998
0.949
1
0.988
Specificity
1
0.976
0.998
0.996
0.807
0.981
k
=
200
k
=
400
# incorrect
64
29
32
210
141
130
Sensitivity
0.682
0.968
0.847
0.479
0.948
0.696
Specificity
1
0.988
0.999
0.999
0.925
0.995
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
100
# incorrect
0
1
1
4
9
2
Sensitivity
0.975
0.994
0.993
0.961
0.998
0.99
Specificity
1
1
1
1
0.995
0.999
k
=
200
k
=
400
# incorrect
15
34
7
82
214
50
Sensitivity
0.928
0.998
0.98
0.803
0.994
0.924
Specificity
1
0.981
0.998
0.998
0.868
0.988
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.Table 6
Simulation results for selectingk important variables among 2,000 candidates including the most and least important surrogate variables across 600 subjects.
Bonferroni
FDR
TT
n.sv = 5
n.sv = 10
n.sv = 15
n.sv = 5
n.sv = 10
n.sv = 15
n.sv = 5
n.sv = 10
n.sv = 15
Most important surrogate variables included
# incorrect
k
=
10
0
0
0
6
6
6
5
4
3
k
=
100
2
3
2
19
19
19
7
5
5
k
=
200
6
6
5
26
29
29
6
6
4
k
=
400
12
11
11
40
38
39
7
7
7
Sensitivity
k
=
10
1
1
1
1
1
1
1
1
1
k
=
100
0.98
0.97
0.98
1
1
1
0.98
0.98
0.98
k
=
200
0.97
0.97
0.975
0.995
0.995
0.995
0.985
0.985
0.985
k
=
400
0.97
0.973
0.973
1
1
1
0.988
0.988
0.985
Specificity
k
=
10
1
1
1
0.997
0.997
0.997
0.997
0.998
0.998
k
=
100
1
1
1
0.99
0.99
0.99
0.997
0.998
0.998
k
=
200
1
1
1
0.986
0.984
0.984
0.998
0.998
0.999
k
=
400
1
1
1
0.975
0.976
0.976
0.999
0.999
0.999
Most important surrogate variables not included
# incorrect
k
=
10
10
10
10
10
10
10
10
10
10
k
=
100
100
100
100
100
100
100
100
100
100
k
=
200
200
200
200
200
200
200
200
200
200
k
=
400
400
400
400
400
400
400
400
400
400
Sensitivity
k
=
10
0
0
0
0
0
0
0
0
0
k
=
100
0
0
0
0
0
0
0
0
0
k
=
200
0
0
0
0
0
0
0
0
0
k
=
400
0
0
0
0
0
0
0
0
0
Specificity
k
=
10
1
1
1
1
1
1
1
1
1
k
=
100
1
1
1
1
1
1
1
1
1
k
=
200
1
1
1
1
1
1
1
1
1
k
=
400
1
1
1
1
1
1
1
1
1
FDR: false discovery rate, TT: training and testing,n.sv = number of surrogate variables, and k: the number of truly important CpG sites out of 2,000.We further examined confounding effect of surrogate variables. To achieve this, we followed the same simulation scenarios noted earlier but, instead of assuming independence between all the (observed and unobserved) independent variables, we allow the 5 unobserved variables to be correlated with observed variablex
1 with a correlation of 0
.
7
i
-
j where i
,
j
∈
L and L
=
(
0,1
,
2,3
,
4,5
) represents an arbitrary locations of the variables. For example, the observed variable is identified at location 0 and the 5 unobserved variables are locations (
1,2
,
3,4
,
5
). Then we generated 100 Monte Carlo replicates based on normal distributions with a sample size of n
=
600. In the simulations, we considered the impact of including and excluding identified surrogate variables. The screening results (results not shown) are consistent with the previous findings when no correlations between unobserved and observed variables were assumed (Tables 1 and 2); that is, including the identified surrogate variable substantially improves the screening statistics (number of incorrectness, sensitivity, and specificity).Based on the simulations, it is advised that users follow default setting for the TT screening method:2
/
3 of the data for training, “two-step” for SVA analysis as described in the published literature [5], 100 iterations for the total number of TT screenings, 50% as the cutoff proportion of those 100 iterations, and 0.05 significance level for the training and testing data. We recommend 100 iterations as a balance between computing efficiency and adequate resampling of the subjects to decipher true associations. We recommend 50% as the default because in general informative CpGs are sparse compared to the number of candidate CpG sites, in which case, as seen in simulations the 50% cutoff percentage is suitable for small and large sample sizes.


### 3.3. Real Data Analysis

We applied thettScreening package to 383,998 CpG sites with DNA methylation available for 245 subjects. The data were collected from a study cohort of 18-year-old females on the Isle of Wight (IOW) in the United Kingdom [17]. We examined a single factor that may be potentially associated with DNA methylation, maternal smoking status during pregnancy (0/1). It is thought that maternal smoking in utero increases chances of asthma and wheezing in children and modifies defensive mechanisms such as xenobiotic detoxification systems, antioxidant responses, and damage repair mechanisms [18]. Certain modifications to such systems have been known to increase risk of lung cancer [18].We applied the TT, FDR-based, and Bonferroni-based methods to identify potentially important CpG sites. In the TT method, 2/3 of the samples were used for training and the remaining for testing. The number of iterations was set at 100 with relative frequency of 50% as the cutoff and the significance level in both training and testing steps was set at 0.05. The significance level for FDR-based and Bonferroni-based methods was set to 0.05. OLS regression was used to estimate associations andp values. Other settings were chosen as the default.FDR- and Bonferroni-based methods, respectively, identified ten and five CpG sites associated with maternal smoking status. The five CpG sites identified by the Bonferroni-based method were also included in the ten CpG sites identified by FDR. The TT screening method identified 91 CpG sites potentially linked to maternal smoking status. The 10 CpG sites collectively identified by Bonferroni and FDR were included in the 91 identified by TT. Significant CpG site locations and annotated genes identified are included in Table7 in Figure 6.Table 7
Information for CpG sites associated with maternal smoking identified by the training and testing screening method.
Model 2: maternal smoking only
CpG site
Chromosome: map
Gene
Promoter region
CpG island
name
information
name
and exon
location
cg20464068
chr1: 14026054
PRDM2
TSS1500
14026481–14027200
cg07951355
chr1: 40123717
cg14179389
chr1: 92947961
GFI1
Body
92945907–92952609
cg00909806
chr1: 212688762
212688031–212688448
cg26143053
chr2: 3718125
ALLC
5′UTR
cg03690080
chr2: 39188006
LOC100271715, LOC375196
Body, TSS1500
39186777–39187968
cg19243656
chr2: 73340094
RAB11FIP5
5′UTR, 1stExon
73339292–73340733
cg18703066‡
chr2: 105363536
cg07516970
chr2: 157181359
NR4A2
3′UTR
157184389–157184632
cg19273101
chr2: 191734576
cg14075934
chr2: 200137014
SATB2
Body
cg03158780
chr3: 577964
cg10663973
chr4: 6642559
MRFAP1
1stExon, 5′UTR
6642194–6643322
cg10364374
chr4: 125857940
cg23743778
chr4: 140357491
140357362–140357610
cg21401642
chr4: 174421114
174421347–174421559
cg17924476†
chr5: 323794
AHRR
Body
320788–323010
cg05575921
chr5: 373378
AHRR
Body
373842–374426
cg21161138
chr5: 399360
AHRR
Body
cg12287936
chr5: 1800606
NDUFS6, MRPL36
TSS1500
1799461–1801905
cg16244648
chr5: 141555043
cg18349863
chr6: 29912713
HLA-A
Body
29910202–29911367
cg11492288
chr6: 30290596
HCG18
Body
30294169–30295071
cg04325960
chr6: 147124986
LOC729176, C6orf103
TSS200, Body
cg00004963
chr6: 147124996
LOC729176, C6orf103
TSS200, Body
cg11881038
chr6: 154408701
OPRM1
Body, 1stExon, 5′UTR
cg20418529
chr6: 166260012
cg00794911
chr6: 166260532
cg18132363†
chr6: 166260572
cg08634229
chr6: 169326603
cg06769202
chr7: 27142535
HOXA2
TSS200
27143181–27143479
cg23206851
chr7: 27143046
HOXA2
TSS1500
27143181–27143479
cg02225599
chr7: 27143252
HOXA2
TSS1500
27143181–27143479
cg10319053
chr7: 27143370
HOXA2
TSS1500
27143181–27143479
cg00445443
chr7: 27143478
HOXA2
TSS1500
27143181–27143479
cg06401979
chr7: 27143717
HOXA2
TSS1500
27143181–27143479
cg11986226
chr7: 40026390
CDK13
Body
cg19089201‡
chr7: 45002287
MYO1G
3′UTR
45002111–45002845
cg04180046‡
chr7: 45002736
MYO1G
Body
45002111–45002845
cg12803068‡
chr7: 45002919
MYO1G
Body
45002111–45002845
cg25949550‡
chr7: 145814306
CNTNAP2
Body
145813030–145814084
cg11207515
chr7: 146904205
CNTNAP2
Body
cg21015808
chr7: 149809179
cg21330896
chr7: 28205705
ZNF395
3′UTR
cg17199018
chr7: 28206278
ZNF395
Body
cg04690729
chr8: 133494328
KCNQ3
TSS1500
133492398–133493586
cg15707110
chr8: 144311102
144311708–144311985
cg08126560
chr9: 92291523
LOC100129066
Body
92291268–92291524
cg13393408
chr9: 132874232
GPR107
Body
cg11813497
chr10: 14372879
FRMD4A
TSS200
cg12490835
chr10: 22623821
22623350–22625875
cg26520012
chr10: 42672589
42672509–42673432
cg05329352
chr10: 112838983
ADRA2A
1stExon
112835990–112839303
cg18424850
chr10: 132945786
TCERG1L
Body
cg19494188
chr11: 1466780
BRSK2
Body
1466304–1467210
cg26204383
chr11: 2435667
TRPM5
Body
2435295–2436651
cg14436038
chr11: 6494706
TRIM3
5′UTR
6494725–6495453
cg15627089
chr11: 16625751
16626053–16629180
cg25160605
chr11: 21087846
NELL1
Body
cg17517598
chr11: 61659090
FADS3
TSS200
61658569–61659592
cg10788371
chr11: 76381040
LRRC32
5′UTR, 1stExon
76381449–76382295
cg11395306
chr11: 98939366
CNTN5
5′UTR
cg01186919
chr11: 111742365
ALG9
TSS1500, TSS200
111741953–111742292
cg02820646
chr11: 115398838
cg05730269
chr11: 118477055
PHLDB1
TSS200, TSS1500
118478235–118481896
cg18493761
chr11: 125386885
cg09932758
chr12: 58022542
B4GALNT1
Body
58021294–58022037
cg09644707
chr12: 114885161
114885105–114885418
cg27103591
chr12: 124809023
NCOR2
3′UTR
124808972–124809176
cg02032696
chr14: 67982198
TMEM229B
TSS200
67981514–67982380
cg21511816
chr14: 76597824
76597648–76597911
cg24874277
chr15: 33211107
FMN1
Body
cg03643241
chr15: 44487910
FRMD5
TSS1500
44486741–44487860
cg20596162
chr15: 45408861
DUOXA2
Body
45408573–45409528
cg16754378
chr15: 57179394
LOC145783
Body
57179277–57179838
cg06899985
chr15: 65689298
IGDCC4
Body
65689142–65689362
cg05549655†
chr15: 75019143
CYP1A1
TSS1500
75018186–75019336
cg17852385
chr15: 75019188
CYP1A1
TSS1500
75018186–75019336
cg11924019†
chr15: 75019283
CYP1A1
TSS1500
75018186–75019336
cg18092474†
chr15: 75019302
CYP1A1
TSS1500
75018186–75019336
cg01060282
chr16: 17033575
cg07675285
chr16: 27121267
27121011–27121241
cg11705699
chr16: 87742845
KLHDC4
Body
87742556–87743109
cg09554007
chr16: 89627174
RPL13, SNORD68
5′UTR, 1stExon, TSS1500
89626644–89627869
cg16483033
chr17: 1090441
ABR
1stExon, 5′UTR
cg08682866
chr17: 13818946
cg17624073
chr17: 79393583
BAHCC1
Body
79393341–79393742
cg13723693
chr17: 80656731
RAB40B
TSS200
80655335–80657183
cg18495341
chr17: 80909836
B3GNTL1
Body
cg18449879
chr19: 16045054
CYP4F11
1stExon, Body
cg02543506
chr20: 33876594
FAM83C
Body
33879904–33880215
CpG: cytosine-phosphate-guanine dinucleotide, chr: chromosome, locations of each CpG are for v37 of the human genome,†CpG sites identified by false discovery rate, and ‡CpG sites identified by false discovery rate and Bonferroni-based methods.Figure 6
Top pathways enriched among genes identified by training and testing method from the Isle of Wight data include aryl hydrocarbon receptor signaling and xenobiotic metabolism signaling, which overlap with other pathways including the nicotine degradation II and nicotine degradation III pathways, as well as pathways for the metabolism of cigarette smoke components and producing known effects of smoke exposure.To understand the biological meaning of the identified sets of CpGs, pathway analysis was used. The genes annotated to each significant CpG site were extracted from the 450 K array manifest file (v1.2; available:http://support.illumina.com/downloads/humanmethylation450_15017482_v1-2_product_files.html). Where a CpG site was annotated to more than one gene, all annotated genes were included. The resulting gene lists were analyzed with Ingenuity Pathway Analysis tool (IPA; Qiagen). Statistical significance for each pathway is reported by IPA using p values. The set of 91 CpGs differentially methylated with maternal smoking identified by TT (includes the 5 and 10 CpG sites identified by Bonferroni and FDR, resp.) were mapped to 54 unique genes.The 91 CpGs differentially methylated with maternal smoking included 18 previously identified top maternal smoking-associated CpG sites inAHRR,CYP1A1,GFI1,MYO1G,CNTNAP2 [4, 19–21],FRMD4A [4, 21],LRRC32, and intergenic CpGs nearLOC284998 andPDE10A-SDIM1 [21]. Bonferroni identified five maternal smoking-associated CpGs and FDR identified ten (including all five identified by Bonferroni), all of which have been previously published in association with maternal smoking. However, the TT method identified many more CpGs that have been previously identified as statistically significant in other cohorts, suggesting that these are not simply false positive results, and represent additional truly maternal smoking associated genes worthy of future investigation and validation. Our observation of differential methylation at age 18 in response to maternal smoking during pregnancy also agrees with previous observations that these epigenetic responses are preserved at least into adolescence [20]. The top pathways enriched among genes containing the 91 CpG sites included aryl hydrocarbon receptor signaling (p
=
0.005) and xenobiotic metabolism signaling (p
=
0.030), which support a long-lasting effect of maternal smoking on metabolic pathways controlling responses to smoke exposure. Furthermore, these pathways overlap with other pathways including the nicotine degradation II and nicotine degradation III pathways, as well as pathways for the metabolism of cigarette smoke components, and producing known effects of smoke exposure (Figure 6), such as melatonin degradation [22]. The TT method also identified six differentially methylated CpG sites inHOXA2, whose mutation causes cleft palate [23], the risk of which is known to be affected by maternal smoking [21].In conclusion, the CpGs detected by TT in association with maternal smoking were enriched among pathways related to the known biology of those processes. Disruption or modification of some of these pathways results in greater risk of lung cancer [18]. TT also identified CpGs located in genes not previously identified which indicates potentially new findings.


## 3.1. Simulations

Simulation Scenarios. In total, 2,000 CpG sites across n
=
600 subjects were simulated. In these 2,000 CpG sites, k sites were assumed to be important. Different settings of k were considered, k
=
10, 100, 200, and 400. Among the k important sites, DNA methylation at 90
% of the k CpG sites was associated with two variables x
1 and x
2, their interaction, and 5 unobservable independent uncorrelated variables, and the remaining 10
% of the k sites were associated with x
1, the interaction of x
1 and x
2, and the 5 unobservable independent uncorrelated variables. Variable x
1 is generated from normal distribution with mean 1 and variance 1, and x
2 is a four-level categorical variable generated from multinomial with parameters n and π
=
{
0.15,0.25,0.25,0.35
}. The remaining CpG sites (2000
-
k) were only associated with the 5 unobservable variables and were deemed as unimportant ones. Linear regressions were applied to simulate the data. To assess the robustness of each method (TT screening and FDR- and Bonferroni-based methods), we considered two types of random error in the regressions, one following normal distribution with mean 0 and variance σ
2
=
1.5 and the other χ
2 distribution with a degree freedom of 1. Combining the choices on k and the distributions of random error, in total, we have 8 settings. For each setting, we generated 100 Monte Carlo (MC) replicates. The results presented below include means of the number of incorrect selections (rounded to the nearest integer), estimates of sensitivity, and estimates of specificity across the 100 MC replicates. The number of incorrect selections refers to the number of CpG sites misidentified (refers to both false positive and false negative CpG sites) out of the 2,000 CpG sites.


## 3.2. Results

Variablesx
1 or x
2 would be selected if their interaction effect was statistically significant. All package options were chosen as the default values with training and testing significance levels both set to 0.05. The screening results (Table 1) indicate that, in general, the sensitivity from the FDR-based method is comparable to that from the TT screening method but its specificity is lower when the number of important variables is not sparse; for example, k
=
200. Compared to the Bonferroni-based method, the TT method in general gave better sensitivity and comparable specificity. These were as expected, as the Bonferroni-based method lacks the ability to control type II error while the FDR-based method cannot control well type I error [16]. The TT screening method, on the other hand, has the potential to control both types I and II errors. This is reflected by the results that the TT screening method overall produced the smallest number of incorrectly identified variables, a statistic incorporating information from both sensitivity and specificity. We also performed the screening using robust regressions and similar results were obtained (not shown). These findings are invariant to the distribution pattern of random errors, normally, or skew distributed.Table 1
Simulation results for selectingk important variables among 2,000 candidates with surrogate variables included.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
100
# incorrect
0
1
1
3
8
2
Sensitivity
0.981
0.991
0.994
0.968
0.999
0.991
Specificity
1
1
1
1
0.996
1
k
=
200
k
=
400
# incorrect
13
32
7
78
213
46
Sensitivity
0.936
0.998
0.982
0.814
0.996
0.931
Specificity
1
0.982
0.998
0.998
0.868
0.989
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
100
# incorrect
1
1
1
16
9
7
Sensitivity
0.892
0.958
0.963
0.837
0.982
0.941
Specificity
1
1
1
1
0.996
1
k
=
200
k
=
400
# incorrect
46
29
21
169
158
95
Sensitivity
0.771
0.981
0.906
0.582
0.969
0.790
Specificity
1
0.986
0.999
0.999
0.909
0.993
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.Recall that the default value of the cutoff percentage required for a CpG site to be treated as an informative site across all iterations was 50%, and the significance levels for both the training and testing data set were at 0.05. To further evaluate how the choices of this cutoff percentage and the significance levels in the training and testing steps influence the screening results, we chose a set of the cutoff percentage values ranging from 30 to 90 and set the training significance level to 0.1 instead of default value 0.05. All other settings were kept the same. As seen in Figure2, across the different number of important CpG sites and sample sizes, overall taking the cutoff percentage close to 50% works the best with significance level of 0.05 for both the training and testing steps. If a higher significance level is chosen in the testing step, then a higher value for cutoff percentage should be used.Figure 2
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites. For the TT screening method, two sets of significance levels are considered: (1) 0.05 for training data and 0.1 for testing data; (2) 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000The above examination on cutoff percentage is with sample sizen
=
600. To assess whether and how sample size influences the choice of cutoff percentage, we repeated the above analyses for different sample sizes (and for each sample size, 100 MC replicates were generated), n
=
200, 400, and 800 with all options inttScreening() set at default values except for cutoff percentage. The results (Figures 3–5) indicated that if the number of important variables is expected to be sparse (e.g., 0.5% of the total number variables), the default cutoff value (50%) is a reasonable choice regardless of the sample size. On the other hand, if important variables are not sparse, for example, at least 5% of the total number of variables, then the cutoff value tends to be influenced by the sample size only if sample size is not large, for example, ≤400. In this case, the smaller the sample size is, the lower the cutoff percentage should be taken. This is as expected; when the sample size is smaller, the probability of true positives will be lower. We thus do not expect a higher proportion of correct selection and consequently a lower cutoff percentage is desired. From our simulations, we recommend 20% cutoff if both the following two conditions are not satisfied: (1) important variables are not expected to be sparse (which rarely happens in high-throughput and high-dimensional data) and (2) we have small samples compared to the number of candidate CpG sites. However, if sample size is large, the default cutoff 50% is still recommended. The selection results for sample sizes n
=
200, 400, and 800 following these recommendations are included in the appendix (Tables 3, 4, and 5).Figure 3
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites across 200 subjects. For the TT screening method, significance levels considered are 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000Figure 4
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites across 400 subjects. For the TT screening method, significance levels considered are 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000Figure 5
Numbers of misidentified CpG sites versus cutoff frequency for a CpG site being potentially important (based on ordinary least squares regressions). The true numbers of important CpG sites are (a) 10, (b) 100, (c) 200, and (d) 400 out of 2,000 CpG sites across 800 subjects. For the TT screening method, significance levels considered are 0.05 for both training and testing data. For the FDR-based and Bonferroni methods, the level was set at 0.05.
(a)
10 important CpG sites out of 2000
(b)
100 important CpG sites out of 2000
(c)
200 important CpG sites out of 2000
(d)
400 important CpG sites out of 2000To demonstrate the benefit of using the identified surrogate variables in screening, we reanalyzed the simulated data (forn
=
600) but with the surrogate variable analysis excluded (this option is available in the package). The relative patterns of statistics between different methods were intact compared to the results when surrogate variables were included in the screening process. However, when surrogate variables were not considered, a drastic drop in the mean sensitivity was observed across all the settings (Table 2). The sensitivity measures across all the 100 replicates ranged from 0 to 30%, implying that all methods had trouble deciphering CpG sites that are truly important. However, instead of completely excluding surrogate variables, what if we include a set of surrogate variables with each surrogate variable explaining quite a small portion of variation (i.e., less informative)? To examine this, we extracted surrogate variables that are less informative and reanalyzed the data. Similar findings as shown in Table 2 were concluded (Table 6), indicating the importance of including informative surrogate variables in order to substantially improve the quality of selection.Table 2
Simulation results for selectingk important variables among 2,000 candidates with surrogate variables excluded.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
100
# incorrect
10
51
9
99
114
90
Sensitivity
0.011
0.054
0.100
0.014
0.167
0.102
Specificity
1
0.979
1
1
0.984
1
k
=
200
k
=
400
# incorrect
197
182
180
394
310
359
Sensitivity
0.015
0.232
0.102
0.015
0.288
0.104
Specificity
1
0.984
1
1
0.984
1
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
100
# incorrect
10
38
9
98
103
91
Sensitivity
0.018
0.053
0.095
0.015
0.158
0.096
Specificity
1
0.986
1
1
0.99
1
k
=
200
k
=
400
# incorrect
197
176
181
394
314
363
Sensitivity
0.014
0.217
0.094
0.014
0.273
0.095
Specificity
1
0.989
1
1
0.986
1
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.Table 3
Simulation results for selectingk important variables among 2,000 candidates including surrogate variables across 200 subjects.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
10
0
†
# incorrect
7
6
6
77
43
41
Sensitivity
0.292
0.423
0.475
0.227
0.605
0.641
Specificity
1
1
1
1
0.998
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
165
87
93
361
214
244
Sensitivity
0.175
0.616
0.571
0.099
0.565
0.432
Specificity
1
0.994
0.996
1
0.975
0.989
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
10
0
†
# incorrect
9
8
7
88
62
55
Sensitivity
0.126
0.193
0.307
0.124
0.408
0.498
Specificity
1
1
1
1
0.998
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
182
123
119
380
274
288
Sensitivity
0.093
0.422
0.435
0.051
0.377
0.314
Specificity
1
0.996
0.997
1
0.985
0.992
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.†Cutoff percentage of 20%.Table 4
Simulation results for selectingk important variables among 2,000 candidates including surrogate variables across 400 subjects.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
10
0
†
# incorrect
2
1
1
24
10
9
Sensitivity
0.851
0.949
0.941
0.761
0.968
0.965
Specificity
1
1
1
1
0.997
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
64
29
21
210
141
93
Sensitivity
0.682
0.968
0.947
0.479
0.948
0.864
Specificity
1
0.988
0.994
0.999
0.925
0.976
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
10
0
†
# incorrect
4
3
3
47
19
18
Sensitivity
0.612
0.767
0.786
0.531
0.878
0.883
Specificity
1
1
1
1
0.997
0.997
k
=
20
0
†
k
=
40
0
†
# incorrect
111
42
42
289
146
141
Sensitivity
0.445
0.88
0.837
0.281
0.841
0.717
Specificity
1
0.99
0.995
0.999
0.949
0.983
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.†Cutoff percentage of 20%.Table 5
Simulation results for selectingk important variables among 2,000 candidates including surrogate variables across 800 subjects.
Statistics
Bon
FDR
TT
Bon
FDR
TT
Random error normally distributed
k
=
10
k
=
100
# incorrect
2
44
5
26
309
36
Sensitivity
0.991
1
0.998
0.949
1
0.988
Specificity
1
0.976
0.998
0.996
0.807
0.981
k
=
200
k
=
400
# incorrect
64
29
32
210
141
130
Sensitivity
0.682
0.968
0.847
0.479
0.948
0.696
Specificity
1
0.988
0.999
0.999
0.925
0.995
Random errorχ
2 distributed (df = 1)
k
=
10
k
=
100
# incorrect
0
1
1
4
9
2
Sensitivity
0.975
0.994
0.993
0.961
0.998
0.99
Specificity
1
1
1
1
0.995
0.999
k
=
200
k
=
400
# incorrect
15
34
7
82
214
50
Sensitivity
0.928
0.998
0.98
0.803
0.994
0.924
Specificity
1
0.981
0.998
0.998
0.868
0.988
Bon: Bonferroni, FDR: false discovery rate, and TT: training and testing.Table 6
Simulation results for selectingk important variables among 2,000 candidates including the most and least important surrogate variables across 600 subjects.
Bonferroni
FDR
TT
n.sv = 5
n.sv = 10
n.sv = 15
n.sv = 5
n.sv = 10
n.sv = 15
n.sv = 5
n.sv = 10
n.sv = 15
Most important surrogate variables included
# incorrect
k
=
10
0
0
0
6
6
6
5
4
3
k
=
100
2
3
2
19
19
19
7
5
5
k
=
200
6
6
5
26
29
29
6
6
4
k
=
400
12
11
11
40
38
39
7
7
7
Sensitivity
k
=
10
1
1
1
1
1
1
1
1
1
k
=
100
0.98
0.97
0.98
1
1
1
0.98
0.98
0.98
k
=
200
0.97
0.97
0.975
0.995
0.995
0.995
0.985
0.985
0.985
k
=
400
0.97
0.973
0.973
1
1
1
0.988
0.988
0.985
Specificity
k
=
10
1
1
1
0.997
0.997
0.997
0.997
0.998
0.998
k
=
100
1
1
1
0.99
0.99
0.99
0.997
0.998
0.998
k
=
200
1
1
1
0.986
0.984
0.984
0.998
0.998
0.999
k
=
400
1
1
1
0.975
0.976
0.976
0.999
0.999
0.999
Most important surrogate variables not included
# incorrect
k
=
10
10
10
10
10
10
10
10
10
10
k
=
100
100
100
100
100
100
100
100
100
100
k
=
200
200
200
200
200
200
200
200
200
200
k
=
400
400
400
400
400
400
400
400
400
400
Sensitivity
k
=
10
0
0
0
0
0
0
0
0
0
k
=
100
0
0
0
0
0
0
0
0
0
k
=
200
0
0
0
0
0
0
0
0
0
k
=
400
0
0
0
0
0
0
0
0
0
Specificity
k
=
10
1
1
1
1
1
1
1
1
1
k
=
100
1
1
1
1
1
1
1
1
1
k
=
200
1
1
1
1
1
1
1
1
1
k
=
400
1
1
1
1
1
1
1
1
1
FDR: false discovery rate, TT: training and testing,n.sv = number of surrogate variables, and k: the number of truly important CpG sites out of 2,000.We further examined confounding effect of surrogate variables. To achieve this, we followed the same simulation scenarios noted earlier but, instead of assuming independence between all the (observed and unobserved) independent variables, we allow the 5 unobserved variables to be correlated with observed variablex
1 with a correlation of 0
.
7
i
-
j where i
,
j
∈
L and L
=
(
0,1
,
2,3
,
4,5
) represents an arbitrary locations of the variables. For example, the observed variable is identified at location 0 and the 5 unobserved variables are locations (
1,2
,
3,4
,
5
). Then we generated 100 Monte Carlo replicates based on normal distributions with a sample size of n
=
600. In the simulations, we considered the impact of including and excluding identified surrogate variables. The screening results (results not shown) are consistent with the previous findings when no correlations between unobserved and observed variables were assumed (Tables 1 and 2); that is, including the identified surrogate variable substantially improves the screening statistics (number of incorrectness, sensitivity, and specificity).Based on the simulations, it is advised that users follow default setting for the TT screening method:2
/
3 of the data for training, “two-step” for SVA analysis as described in the published literature [5], 100 iterations for the total number of TT screenings, 50% as the cutoff proportion of those 100 iterations, and 0.05 significance level for the training and testing data. We recommend 100 iterations as a balance between computing efficiency and adequate resampling of the subjects to decipher true associations. We recommend 50% as the default because in general informative CpGs are sparse compared to the number of candidate CpG sites, in which case, as seen in simulations the 50% cutoff percentage is suitable for small and large sample sizes.


## 3.3. Real Data Analysis

We applied thettScreening package to 383,998 CpG sites with DNA methylation available for 245 subjects. The data were collected from a study cohort of 18-year-old females on the Isle of Wight (IOW) in the United Kingdom [17]. We examined a single factor that may be potentially associated with DNA methylation, maternal smoking status during pregnancy (0/1). It is thought that maternal smoking in utero increases chances of asthma and wheezing in children and modifies defensive mechanisms such as xenobiotic detoxification systems, antioxidant responses, and damage repair mechanisms [18]. Certain modifications to such systems have been known to increase risk of lung cancer [18].We applied the TT, FDR-based, and Bonferroni-based methods to identify potentially important CpG sites. In the TT method, 2/3 of the samples were used for training and the remaining for testing. The number of iterations was set at 100 with relative frequency of 50% as the cutoff and the significance level in both training and testing steps was set at 0.05. The significance level for FDR-based and Bonferroni-based methods was set to 0.05. OLS regression was used to estimate associations andp values. Other settings were chosen as the default.FDR- and Bonferroni-based methods, respectively, identified ten and five CpG sites associated with maternal smoking status. The five CpG sites identified by the Bonferroni-based method were also included in the ten CpG sites identified by FDR. The TT screening method identified 91 CpG sites potentially linked to maternal smoking status. The 10 CpG sites collectively identified by Bonferroni and FDR were included in the 91 identified by TT. Significant CpG site locations and annotated genes identified are included in Table7 in Figure 6.Table 7
Information for CpG sites associated with maternal smoking identified by the training and testing screening method.
Model 2: maternal smoking only
CpG site
Chromosome: map
Gene
Promoter region
CpG island
name
information
name
and exon
location
cg20464068
chr1: 14026054
PRDM2
TSS1500
14026481–14027200
cg07951355
chr1: 40123717
cg14179389
chr1: 92947961
GFI1
Body
92945907–92952609
cg00909806
chr1: 212688762
212688031–212688448
cg26143053
chr2: 3718125
ALLC
5′UTR
cg03690080
chr2: 39188006
LOC100271715, LOC375196
Body, TSS1500
39186777–39187968
cg19243656
chr2: 73340094
RAB11FIP5
5′UTR, 1stExon
73339292–73340733
cg18703066‡
chr2: 105363536
cg07516970
chr2: 157181359
NR4A2
3′UTR
157184389–157184632
cg19273101
chr2: 191734576
cg14075934
chr2: 200137014
SATB2
Body
cg03158780
chr3: 577964
cg10663973
chr4: 6642559
MRFAP1
1stExon, 5′UTR
6642194–6643322
cg10364374
chr4: 125857940
cg23743778
chr4: 140357491
140357362–140357610
cg21401642
chr4: 174421114
174421347–174421559
cg17924476†
chr5: 323794
AHRR
Body
320788–323010
cg05575921
chr5: 373378
AHRR
Body
373842–374426
cg21161138
chr5: 399360
AHRR
Body
cg12287936
chr5: 1800606
NDUFS6, MRPL36
TSS1500
1799461–1801905
cg16244648
chr5: 141555043
cg18349863
chr6: 29912713
HLA-A
Body
29910202–29911367
cg11492288
chr6: 30290596
HCG18
Body
30294169–30295071
cg04325960
chr6: 147124986
LOC729176, C6orf103
TSS200, Body
cg00004963
chr6: 147124996
LOC729176, C6orf103
TSS200, Body
cg11881038
chr6: 154408701
OPRM1
Body, 1stExon, 5′UTR
cg20418529
chr6: 166260012
cg00794911
chr6: 166260532
cg18132363†
chr6: 166260572
cg08634229
chr6: 169326603
cg06769202
chr7: 27142535
HOXA2
TSS200
27143181–27143479
cg23206851
chr7: 27143046
HOXA2
TSS1500
27143181–27143479
cg02225599
chr7: 27143252
HOXA2
TSS1500
27143181–27143479
cg10319053
chr7: 27143370
HOXA2
TSS1500
27143181–27143479
cg00445443
chr7: 27143478
HOXA2
TSS1500
27143181–27143479
cg06401979
chr7: 27143717
HOXA2
TSS1500
27143181–27143479
cg11986226
chr7: 40026390
CDK13
Body
cg19089201‡
chr7: 45002287
MYO1G
3′UTR
45002111–45002845
cg04180046‡
chr7: 45002736
MYO1G
Body
45002111–45002845
cg12803068‡
chr7: 45002919
MYO1G
Body
45002111–45002845
cg25949550‡
chr7: 145814306
CNTNAP2
Body
145813030–145814084
cg11207515
chr7: 146904205
CNTNAP2
Body
cg21015808
chr7: 149809179
cg21330896
chr7: 28205705
ZNF395
3′UTR
cg17199018
chr7: 28206278
ZNF395
Body
cg04690729
chr8: 133494328
KCNQ3
TSS1500
133492398–133493586
cg15707110
chr8: 144311102
144311708–144311985
cg08126560
chr9: 92291523
LOC100129066
Body
92291268–92291524
cg13393408
chr9: 132874232
GPR107
Body
cg11813497
chr10: 14372879
FRMD4A
TSS200
cg12490835
chr10: 22623821
22623350–22625875
cg26520012
chr10: 42672589
42672509–42673432
cg05329352
chr10: 112838983
ADRA2A
1stExon
112835990–112839303
cg18424850
chr10: 132945786
TCERG1L
Body
cg19494188
chr11: 1466780
BRSK2
Body
1466304–1467210
cg26204383
chr11: 2435667
TRPM5
Body
2435295–2436651
cg14436038
chr11: 6494706
TRIM3
5′UTR
6494725–6495453
cg15627089
chr11: 16625751
16626053–16629180
cg25160605
chr11: 21087846
NELL1
Body
cg17517598
chr11: 61659090
FADS3
TSS200
61658569–61659592
cg10788371
chr11: 76381040
LRRC32
5′UTR, 1stExon
76381449–76382295
cg11395306
chr11: 98939366
CNTN5
5′UTR
cg01186919
chr11: 111742365
ALG9
TSS1500, TSS200
111741953–111742292
cg02820646
chr11: 115398838
cg05730269
chr11: 118477055
PHLDB1
TSS200, TSS1500
118478235–118481896
cg18493761
chr11: 125386885
cg09932758
chr12: 58022542
B4GALNT1
Body
58021294–58022037
cg09644707
chr12: 114885161
114885105–114885418
cg27103591
chr12: 124809023
NCOR2
3′UTR
124808972–124809176
cg02032696
chr14: 67982198
TMEM229B
TSS200
67981514–67982380
cg21511816
chr14: 76597824
76597648–76597911
cg24874277
chr15: 33211107
FMN1
Body
cg03643241
chr15: 44487910
FRMD5
TSS1500
44486741–44487860
cg20596162
chr15: 45408861
DUOXA2
Body
45408573–45409528
cg16754378
chr15: 57179394
LOC145783
Body
57179277–57179838
cg06899985
chr15: 65689298
IGDCC4
Body
65689142–65689362
cg05549655†
chr15: 75019143
CYP1A1
TSS1500
75018186–75019336
cg17852385
chr15: 75019188
CYP1A1
TSS1500
75018186–75019336
cg11924019†
chr15: 75019283
CYP1A1
TSS1500
75018186–75019336
cg18092474†
chr15: 75019302
CYP1A1
TSS1500
75018186–75019336
cg01060282
chr16: 17033575
cg07675285
chr16: 27121267
27121011–27121241
cg11705699
chr16: 87742845
KLHDC4
Body
87742556–87743109
cg09554007
chr16: 89627174
RPL13, SNORD68
5′UTR, 1stExon, TSS1500
89626644–89627869
cg16483033
chr17: 1090441
ABR
1stExon, 5′UTR
cg08682866
chr17: 13818946
cg17624073
chr17: 79393583
BAHCC1
Body
79393341–79393742
cg13723693
chr17: 80656731
RAB40B
TSS200
80655335–80657183
cg18495341
chr17: 80909836
B3GNTL1
Body
cg18449879
chr19: 16045054
CYP4F11
1stExon, Body
cg02543506
chr20: 33876594
FAM83C
Body
33879904–33880215
CpG: cytosine-phosphate-guanine dinucleotide, chr: chromosome, locations of each CpG are for v37 of the human genome,†CpG sites identified by false discovery rate, and ‡CpG sites identified by false discovery rate and Bonferroni-based methods.Figure 6
Top pathways enriched among genes identified by training and testing method from the Isle of Wight data include aryl hydrocarbon receptor signaling and xenobiotic metabolism signaling, which overlap with other pathways including the nicotine degradation II and nicotine degradation III pathways, as well as pathways for the metabolism of cigarette smoke components and producing known effects of smoke exposure.To understand the biological meaning of the identified sets of CpGs, pathway analysis was used. The genes annotated to each significant CpG site were extracted from the 450 K array manifest file (v1.2; available:http://support.illumina.com/downloads/humanmethylation450_15017482_v1-2_product_files.html). Where a CpG site was annotated to more than one gene, all annotated genes were included. The resulting gene lists were analyzed with Ingenuity Pathway Analysis tool (IPA; Qiagen). Statistical significance for each pathway is reported by IPA using p values. The set of 91 CpGs differentially methylated with maternal smoking identified by TT (includes the 5 and 10 CpG sites identified by Bonferroni and FDR, resp.) were mapped to 54 unique genes.The 91 CpGs differentially methylated with maternal smoking included 18 previously identified top maternal smoking-associated CpG sites inAHRR,CYP1A1,GFI1,MYO1G,CNTNAP2 [4, 19–21],FRMD4A [4, 21],LRRC32, and intergenic CpGs nearLOC284998 andPDE10A-SDIM1 [21]. Bonferroni identified five maternal smoking-associated CpGs and FDR identified ten (including all five identified by Bonferroni), all of which have been previously published in association with maternal smoking. However, the TT method identified many more CpGs that have been previously identified as statistically significant in other cohorts, suggesting that these are not simply false positive results, and represent additional truly maternal smoking associated genes worthy of future investigation and validation. Our observation of differential methylation at age 18 in response to maternal smoking during pregnancy also agrees with previous observations that these epigenetic responses are preserved at least into adolescence [20]. The top pathways enriched among genes containing the 91 CpG sites included aryl hydrocarbon receptor signaling (p
=
0.005) and xenobiotic metabolism signaling (p
=
0.030), which support a long-lasting effect of maternal smoking on metabolic pathways controlling responses to smoke exposure. Furthermore, these pathways overlap with other pathways including the nicotine degradation II and nicotine degradation III pathways, as well as pathways for the metabolism of cigarette smoke components, and producing known effects of smoke exposure (Figure 6), such as melatonin degradation [22]. The TT method also identified six differentially methylated CpG sites inHOXA2, whose mutation causes cleft palate [23], the risk of which is known to be affected by maternal smoking [21].In conclusion, the CpGs detected by TT in association with maternal smoking were enriched among pathways related to the known biology of those processes. Disruption or modification of some of these pathways results in greater risk of lung cancer [18]. TT also identified CpGs located in genes not previously identified which indicates potentially new findings.


## 4. Conclusions

We developed a unique screening procedure built into an R package for the purpose of screening important variables. It includes the developed method that involves training and testing steps (the TT screening method), along with another two existing methods, the method controlling for FDR and the other controlling overall significance level through the Bonferroni method. Simulations are used to demonstrate and assess the TT screening method.Overall, the TT screening method produced comparable sensitivity results to that of FDR-based method and comparable specificity results of Bonferroni-based methods. However, the number of misidentified CpG sites from the TT method is in general smaller than those from the other two approaches. The findings on sensitivity and specificity were as expected, because the Bonferroni-based method lacks the ability to control type II error while the FDR-based method cannot control well type I error [16]. The TT screening method, on the other hand, has the potential to control both the types I and II errors, which was supported by the smaller numbers of incorrect detections. It was also noticed that, by incorporating surrogate variable analysis, all three methods produced higher sensitivity measures. In the real data application, the TT method identified a larger number of CpGs compared to the other two methods. The subsequent pathway analyses further support the practical strength of the proposed method. It is worth noting that since the screening process is to identify the CpG sites, to properly assess the functionality of these sites, it is critical to evaluate them jointly, for example, using pathway analyses.An important contribution of the developed process is its computing efficiency. Although we combined three methods into one package, the computing is not burdensome. The computing time for screening 383,998 CpG sites (the real data analysis) with thettScreening() function using default options only takes 82 minutes on a personal computer with a central processing unit (CPU) of 2.0 gigahertz (GHz) and memory of 4.0 gigabytes (GB) of random-access memory (RAM). ThettScreening() function always automatically adjusts for multiple testing using three methods, FDR, Bonferroni, and the training and testing method.Lastly, the versatility of the proposed screening method allows it to be applied to a variety of scientific fields in which large data or high dimensionality is a computational problem. DNA methylation is of growing interest in all aspects of public health, in cancer and genetics/epigenetics specifically. TT reduces dimensionality in a timely manner while controlling for types I and II errors and adjusting unknown latent variables estimated using the surrogate variable analysis. Finally, the package incorporates a variety of options which allows the user to create very specific settings while maintaining convenient usability.



---

*Source: 2615348-2016-03-13.xml*