For real data analysis, we chose a Korean GWAS data set collected since 2007 by The Korean Association REsource (KARE) project (Cho et al., 2009). In this project, all participants were recruited from either of two region-based cohorts (rural Ansung and urban Ansan). The total number of participants was 10,038 (5,018 from Ansung and 5,020 from Ansan), and they were all genotyped, using genomic DNA from peripheral blood, using the Affymetrix (Santa Clara, CA, United States) Genome-Wide Human SNP array 5.0, containing 500,568 SNPs. For quality control, we followed the same process used in a previous study (Oh et al., 2016). As a result, we finally obtained 8,842 individuals and 327,872 SNPs, and the processed data set was used in our real data analysis. The study was reviewed and approved by the Institutional Review Board of Seoul National University (IRB No. E1908/001-004).

Suppose that there are N samples, each with a dependent variable Y, and J features X₁,…,X_J, representing features from a multi-omics data set. In general, for a significance test of association between a specific X_j and Y, the null distribution of the test statistic S consists of test statistics from permuted data sets, and we call the statistics s_r, where r = 1,2,…,R, with R denoting the total number of permutation rounds for the feature. Then, the observed value, s_obs (i.e., the original value of the test statistic, S) is compared to the null distribution of S, and the significance is assessed by the proportion of s_r values more extreme than s_obs. For exact generation of the null distribution, N! iterations are required. However, when N! is too large, R iterations of random shuffling (R) (R≪N!) are generally used for assessing computational feasibility in terms of Monte-Carlo estimation. A finding that a s_obs value is larger than the simulated s_r values implies that the test is more supportive of the alternative hypothesis, and the p-value is then calculated by the following equation:

where I(⋅) is an indicator function, and +1 in the numerator and denominator can be omitted.

When the number of features is multiple, the p-value threshold should be adjusted for a multiple testing comparison. For example, a typical p-value threshold is 0.05, and, if there are 1,000 features for association tests, then the p-value threshold becomes 0.05/1,000, for the Bonferroni correction. In other words, when a feature has a p-value smaller than this adjusted p-value threshold it is reported as significant. Therefore, the possibility of I(⋅) = 1 (more extreme than the observed statistic value) is extremely low for this feature. On the other hand, if I(⋅) = 1 frequently appears in a feature, the p-value of the feature may be closer to 1, meaning that it may not be significant and would therefore be of no interest to researchers. Let p_raw be an unadjusted p-value threshold (e.g., 0.05) and p_adj be an adjusted p-value threshold, for each feature, after the multiple testing correction (e.g., 0.05/J by Bonferroni correction). p_adj is then the significance level for which we need to detect significant features, and the decision of whether or not to prune a feature, in any specific round, is based on the hypothesis that:

where p implies the true p-value from the permutation approach. In the hypothesis, the significance level for the test needs to be determined, and we call the threshold p_prun. For the hypothesis test, a binomial test can be used, and, based on p_adj and p_prun, we can set an integer C_prun that satisfies p_prun in a permutation round. Therefore, C_prun is a variable that depends on permutation numbers, while p_adj and p_prun are fixed values for the whole pruning process. Consequently, using this rule, EPNN counts in how many cases a feature has a more extreme test statistic than its observed test statistic value in each permutation round. If a feature is equal to or greater than C_prun in a round, it is removed from the next permutation round. The following is a detailed explanation of the parameter determination.

Let us assume that p_adj = 5 × 10^–5, which is equivalent to a threshold Bonferroni correction with 1,000 features, and p_prun = p_adj. In addition, if we let p_k—r denote a probability of observing at least a number k of test statistics values more extreme than the observed test statistics at the rth permutation round, then pk|r=∑t=kt=r(rt)⁢pa⁢d⁢jt⁢(1-pa⁢d⁢j)r-t. Therefore, if the p-value of a feature is significant, then p_k—r should be equal to or smaller than p_prun. As an illustration, consider the first permutation round. Based on a setting of p_adj = 5 × 10^–5, two probabilities, p_0|1,p_1|1, are given. Because we set p_{prun =}p_adj,p_0|1 will be 1 and p_1—1 will be p_adj, implying that C_prun = 1 is in the first round. For the second round, there are three probabilities, p_0|2,p_1|2 and p_2|2, that can be easily computed. In this case, p_1|2 = 1×10⁻⁴ > p_prunp_2|2 = 10⁻⁹ < p_prun. Therefore C_prun will be 2 for the second round. In this manner, we can obtain C_prun for all permutation rounds conducted. We will show the properties of the parameters in the next section.

In this section, we evaluated the advantages of ENPP compared to a strict permutation approach, including its need for only very few counts for rejecting and removing non-significant features. As a consequence of this attribute, ENPP can greatly reduce total computation time to a feasible level compared to an unpruned permutation approach. To show the desired properties, we artificially generated data sets whose features did not associate with a feature. When the Bonferroni threshold was applied and p_raw = 0.05, the first example had pa⁢d⁢j1 = 0.05/1,000 and the second example had pa⁢d⁢j2 = 0.05/(5 × 10⁵). In addition, we also assumed that p_prun = p_adj for both examples.

Firstly, we investigated the distribution of C_prun values according to each permutation round for pa⁢d⁢j1, and pa⁢d⁢j2, respectively. Using the formula described in the methods, C_prun values were calculated for r = 1,2,…, 10,000, and the resulting values are shown in Figure 1A, which also shows that the values of C_prun for pa⁢d⁢j1 are at most 6 in the 10,000th round. This implies that the threshold is not hard to satisfy and that we can reduce a large proportion of the number of features at each permutation round. In the case of pa⁢d⁢j2, C_prun becomes smaller (Figure 1B). In detail, C_prun is 1 for i = 1, 2 for i ∈ [2, 4,473], and 3 for i ∈ [4,474, 10,000], implying that smaller p_adj values provide smaller C_prun values, although p_prun is proportional to p_adj.

Based on the C_prun values calculated above, we also evaluated the pruned proportion of the total features for each permutation round. Suppose that the p-value of a feature has a uniform distribution, meaning that the feature has no association with a phenotype. In this setting, the pruned proportion of features depends only on C_prun. For example, at the first round, for C_prun(1) = 1, the proportion of pruned features will be ∫01p⁢dp=12. At the second round, for C_prun(2) = 2, no pruning will happen, because the event that C_prun(1) = 1 includes the event that C_prun(2) = 2. At the third round of permutation, for C_prun(3) = 2, the expected pruning proportion after the permutation will be:

In other words, at the first permutation, 12 of the features are expected to be pruned, and 112 of the features are additionally pruned after the third round. In this manner, the expected proportions of remaining features after pruning from 1 to 10,000 permutation rounds are calculated using the C_prun values (Figure 1), and the results are described in Figure 2. Because the cumulative pruning proportion is not easily derived by numerical calculation, we estimated the proportion by simulation using variables from a Bernoulli distribution, with the probability for success taken from a uniform distribution U(0,1). In Figure 2A, only about 2% of features remain after the 100th permutation round in both p_prun settings, thus greatly reducing the number of tests for the data set at the round. However, as C_prun becomes different, the remaining proportions also become different. For example, at the 1000th permutation round, 0.3% of total features remained for pa⁢d⁢j1 and 0.2% for pa⁢d⁢j2. The ratio between the two proportions became larger at the 10,000th permutation round, with 0.057% for the former, pa⁢d⁢j1, and 0.028% for the latter, pa⁢d⁢j2. These results reflect the differences of C_prun provided in Figure 1.

We next assessed computational efficiency by comparing the total permutation time for ENPP to that for the original, unpruned permutation test. The efficiency is represented as a ratio between the number of tests in the original unpruned permutation approach and the cumulative number of tests in the ENPP approach. The total permutation time for a given permutation round in ENPP is calculated by accumulating all permutation times of earlier permutation rounds. Therefore, larger computational efficiencies imply a large timesaving advantage for ENPP analysis. For example, during the first round, there is no reduction of permutation time, but for the second and third permutation rounds, ENPP needs only 12 the computations compared to the original unpruned permutation tests, and 512 the permutations are needed for the fourth round. Therefore, computational efficiency will be 11=1 for the first permutation round, and 1+11+12=43,1+1+11+12+12=32,1+1+1+11+12+12+512=4829. for the second, third, and fourth permutation rounds, respectively. The Inverse Computational Efficiency (ICE) for each permutation round is summarized in Figure 2B. In Figure 2B, ICE does not seem to decrease as fast as the remaining proportion, as shown in Figure 2A, due to the fact that permutation times of precedent rounds accumulate in estimating computational efficiency. Compared to the ordinary unpruned permutation test, only about 7.4% of the computation time is needed at the 100th permutation round in both settings, because they have the same numbers for C_prun and the same resulting remaining proportions. However, as in the remaining proportion of features, ICE became more different in terms of ratios between the two settings as the permutation round progresses. For example, at the 1000th permutation round, ICE is 1.3% for p_prun = 5 × 10^–5 and 1.2% for p_prun = 1 × 10^–7. However, in the 10,000 iteration, 0.23% is needed for the former, p_prun while 0.17% is needed for the latter, p_prun. Thus, the overall computational efficiency improves as the iteration round progresses because the remaining rate of the features grows smaller, and smaller p_prun requires less computation.

On the other hand, we assessed the type 1 error rate of non-associated features from the ENPP approach. For pa⁢d⁢j1 and pa⁢d⁢j2, we generated 10⁶ and 5 × 10⁸ non-associated features from the Bernoulli distribution so that the expected numbers of features with type 1 error are 50 in both settings. We first applied the pruning process to the non-associated features and then the full permutation approach to the remaining unpruned features. After the full permutation approach had been applied, we counted how many non-associated features were found significant at the given significance levels. The type 1 error rates are summarized in Figure 2C, showing that the ENPP approach controls the type 1 error well.

We next applied our approach to a real genome-wide data set (Korea Association REsource: KARE), which has 327,872 SNPs from each of 8,842 individuals (Cho et al., 2009). In order to detect significant SNP features at the Bonferroni significance level in the data set, the ordinary permutation approach (without ENPP) requires at least (1/0.05) × 327,872² = 2.15 × 10¹², a computationally impractical number of tests. Therefore, using a pruning approach for this data set becomes inevitable when the permutation approach is used. For the application of ENPP, we set p_raw = 0.05 and p_prun = p_adj = 0.05/327,872 = 1.52 × 10^–7, and the corresponding C_prun is calculated and described in Figure 3A. Here, we set the number of iterations to 100,000 because simulation analysis found that the remaining proportion of features was 3.7 × 10^–5 at the 100,000th round and the corresponding expected count of remaining features was 3.7 × 10^–5 × 327,872 = 12.13 if all features were assumed not to associate with a phenotype. We selected fasting plasma glucose (FPG) as a phenotype because its distribution is very highly skewed (skewness = 5.32) and the skewness is still high (=2.71) (Kim, 2013) even after log-transformation. Consequently, we expected that this property may produce results that differ between a parametric approach and a permutation approach. For the association analysis, we used age, gender, and living regions as covariates, and we assumed that the genotype of the SNP features has an additive effect on the phenotype. As a test statistic for the permutation test, we used a t-statistic for the genotype effect.

Based on the expected remaining proportion of the features, we found ICE to be 2.4 × 10^–4 at the 100,000th permutation round (Figure 3C), meaning that we needed only 24 times more computation compared to the parametric linear regression approach. This number of permutation tests can be done in a few days, even in a single thread. After implementing the 100,000th iteration of ENPP with the real data set, we plotted the number of remaining features (Figure 3B) and the ICE (Figure 3C) in each round. Those results showed that 46 SNP features remained and that the computational efficiency was 3.7 × 10^–4, implying that some SNP features were candidates for significant features. For each of 46 SNP features, we implemented a 3×10⁷−1 permutation test to provide a p-value not only for Bonferroni correction but also for a genome-wide significance of 5 × 10^–8 (Xu et al., 2014). After implementation of the test, we found that five SNP features passed the Bonferroni threshold, and two SNPs also passed for genome-wide significance (Table 1). On the other hand, the parametric approach found four SNPs for Bonferroni correction, and two SNPs passed genome-wide significance. However, only three SNPs overlapped for the former threshold, and one SNP overlapped for the latter one. To determine substantial differences of p-values between the two approaches, we used an exact binomial test (Clopper and Pearson, 1934) that regarded p-values from the parametric approach as a null hypothesis p-value for the permutation results. From the test, we found that only one SNP (rs7197218G in chromosome 16) showed a significant difference between the two results (Table 1). This SNP showed a more conservative result from the permutation approach; this result may come from type 1 error inflation in the parametric test in the presence of very low minor allele frequency and large differences of variance between FPG values with and without the minor allele (Zimmerman, 2004).