Detecting PCOS susceptibility loci from genome-wide association studies via iterative trend correlation based feature screening.

BACKGROUND: Feature screening plays a critical role in handling ultrahigh dimensional data analyses when the number of features exponentially exceeds the number of observations. It is increasingly common in biomedical research to have case-control (binary) response and an extremely large-scale categorical features. However, the approach considering such data types is limited in extant literature. In this article, we propose a new feature screening approach based on the iterative trend correlation (ITC-SIS, for short) to detect important susceptibility loci that are associated with the polycystic ovary syndrome (PCOS) affection status by screening 731,442 SNP features that were collected from the genome-wide association studies.
RESULTS: We prove that the trend correlation based screening approach satisfies the theoretical strong screening consistency property under a set of reasonable conditions, which provides an appealing theoretical support for its outperformance. We demonstrate that the finite sample performance of ITC-SIS is accurate and fast through various simulation designs.
CONCLUSION: ITC-SIS serves as a good alternative method to detect disease susceptibility loci for clinic genomic data.

Background

Ultrahigh dimensional data with binary response and categorical features has become increasingly prevalent in various fields. Applications using such data exist in genome-wide association studies (GWAS), medical imaging, finance, text mining, among others [1, 2]. The most prevailing gene selection approaches used in genome-wide association studies consider an association of each of the genetic variant using univariate models (i.e, single-SNP models); however, they evaluate the association of each SNP in isolation from the others and hence ignore combined joint effects of multi-loci [3-7]. As a matter of fact, most complex diseases are reported to be mediated through multiple genetic variants, each conferring a small or moderate effect with low penetrance, which obscures the individual significance of each variant [8, 9]. Furthermore, 70−−80% of genomes showing regions of high linkage disequilibrium (LD), which is the nonrandom association of alleles at nearby loci [10-12]. Malo et al. (2008) claimed that single-SNP approaches failed to differentiate truly influential SNPs from spurious SNPs that were merely in LD with the influential SNPs [13]. Therefore, although widely used in GWAS data analyses for its simplicity, single-SNP models have limited power and yield both high false-positive and false-negative results [13-15].

Many joint models can be applied to analyze the association between categorical features and binary response if the dimension is moderate. For example, random forests [16, 17], k-nearest neighbors [18], and support vector machines [19, 20], etc. However, these methods may lose power or become increasingly unstable, and hence intractable, as the dimension of feature space becomes ultrahigh [21]. To fully address the joint effect of multi-loci together with confounding caused by LD for high dimensional data, the penalized regression approaches have been well established and widely used in gene selections [13, 22, 23]. Through intensively investigating 48 different settings with varied LD strength, minor allele frequency, and dimensionality, Michelle et al. empirically verified that the running time and selection success rates decreased dramatically as the dimension of SNPs is greater than 1000 and even worse for 10,000 when applying a multiple logistic penalized ridge regression to gene selection [15]. Furthermore, the computational expediency, statistical accuracy and theoretical disadvantages of penalized regression approaches were concerned for ultrahigh dimensional data analyses [21, 24], which represent the true need for genome-wide association studies. As commented by Fan et al. [21, 25], the fundamental challenges of big data come from the accumulation of aggregate error rates due to a preponderance of noise features. Actually the majority of the information in ultrahigh dimensional data is represented by only a small amount of truly influential features.

Feature screening has garnered considerable attention in recent statistic literature. It filters out a substantial amount of noise features to truly reflect the sparsity principle of the ultrahigh dimensional data. In their seminal paper, Fan and Lv established the underpinnings for what they termed sure independent screening (SIS) and introduced the conceptual framework of the bulk of feature screening literature that come thereafter [26]. Even though many approaches stemming from [26] relaxed the model specification assumptions (see an overview in [27]), many existing SIS-based procedures still tacitly require that the feature variables and the response are continuous [28-31]. Notably, this implicit presupposition of continuity of the variables can be limiting in several application directions, for example in GWAS area.

Motivated by a polycystic ovary syndrome (PCOS) data with 4,099 observations (1,043 cases and 3,056 controls) and 731,442 single nucleotide polymorphisms (SNPs) (p>>n), in this article we propose a new feature screening method based on the iterative trend correlation (we call it ITC-SIS for short). We prove that trend correlation based strong independence screening (TC-SIS) satisfies the strong screening consistency property under a set of reasonable conditions, which is a much stricter criterion than the generally proved sure screening property. Since TC-SIS is a marginal approach, the iterative process on TC-SIS is applied to detect the multi-loci effect each having weak main effects, and separate the individual variants that are truly influential from those confounding spurious variants that are irrelevant to the response but highly correlated with the causative loci due to LD.

There are three existing methods in the statistic feature screening literature that also admit a binary response and categorical features for ultrahigh dimensionality: the maximum marginal likelihood estimator based approach (MMLE-SIS) [32], the distance correlation based approach (DC-SIS) [24], and the Pearson’s chi-squared test based approach (PC-SIS) [2]. We compare TC-SIS with these three most relevant methods and demonstrate that TC-SIS has agreeable accuracy and speed in handling the motivated setting with categorical feature and binary outcome through various finite sample simulation studies. We also demonstrate the ITC-SIS indeed improves TC-SIS through iterative process.

Methods

2.1 Some preliminaries

Trend correlation was applied to measure association between two categorical variables [33, 34]. Specifically, let Y be the response variable and X, j=1,…,p, be the j categorical feature variable. Define be the numeric score assigned to each level of X, where k=1,…,K. Here p is the total number of features and K is the number of levels for X, for which we allow for various number of levels for different features. Let m=0,1 be the encoding of the case-control response Y.

The trend correlation ϱ between the response Y and feature X is defined as follows [35]

where are the frequencies of each level of the feature, response, and individual cell of the contingency table, respectively. We are motivated to utilize ϱ as a feature screening procedure because it possesses the rather salient property of being equal to zero if and only if the two involved variables are independent [35].

2.2 A new independence ranking and screening procedure

Here we describe the details of the proposed TC-SIS screening procedure. For a sample of n observations, we will denote the sample mean score of X by and sample mean of Y as . We then estimate the trend correlation between X and Y as [35]

where , , and represent the sample proportions that are used to estimate the corresponding population proportions , , and p, respectively.

When the features are ordinal, we can interpret as estimating the linear trend between X and Y, e.g., an increase in the observed level of X tends to be associated with decreasing or increasing levels of Y [35]. Therefore, it is suggested that the ordering of and the distance between the scores conform to those of the categorical levels.

The sparsity principle of ultrahigh dimensional data indicates that only a small number of the features truly influence the response. Define as the set of the full model, i.e., all features in the candidate pool. Let be a subset of , i.e., an arbitrary model under consideration.

We use as a marginal utility to rank the importance of each X according to its associations with the response, where higher values correspond to stronger association. Note that is non-negative because the absolute values are used in the numerator. As a output of the TC-SIS feature screening procedure, the selected model are given by

where c is a pre-specified threshold value.

The aim of feature screening is to select the true model or at least select a model that contains the true model. As a matter of further notation, we will denote the true model by and the selected model output from TC-SIS by .

2.3 Theoretical properties

In this section the theoretical properties of the proposed independence screening procedure TC-SIS will be studied. We first define two conditions to facilitate the technical proofs:

Bounds on the standard deviations. Assume that there exists a positive constant σmin such that for all j, This excludes unusual or unreasonable features that are constant and hence have a standard deviation of zero. It should further be noted that an upper bound on σ and σ can also be obtained, by use of Popoviciu’s inequality on variances (see [36]):

Lower bound on the covariance. Assume that ϱ=0 for any . Assume that there exists a positive constant ωmin such that which indicates that the correlation between each truly influential feature and the response is not trivial.

When these two conditions are satisfied, we can establish the following theorems that support the strong screening property for the TC-SIS procedure.

Theorem 1

(Sure Screening Property). Under condition (C1) and removing from (C2) only the assumption that ϱ=0 for any , there exists a positive constant c>0 such that (However, may not converge to 1 as n approaches infinity).

Theorem 2

(Strong Screening Consistency). Given conditions (C1) and (C2), there exists a positive constant c>0 such that

The property of strong screening consistency is much harder to achieve than the (weak) sure screening property because it not only guarantees that the true model is contained in the selected subset, but also ensures that the selected subset is the minimum one containing the true model asymptotically. The proofs of these two theorems are presented in the Supplement file. In addition to the aforementioned two theorems, we also draw two corollaries, which are not themselves related to sure screening, but they are nevertheless important conclusions related to the screening criterion of the TC-SIS method.

Corollary 1

There exists a positive value ϱmin>0 such that for any , we have ϱ>ϱmin>0. This will be shown in Step 1 of the proofs of Theorems 1 and 2.

Corollary 2

The estimator converges uniformly in probability to ϱ. In other words, for any ε>0. This will be shown in Step 2 of the proofs of Theorems 1 and 2.

2.4 Iterative process

Although the proposed TC-SIS approach is powerful at filtering out noise and selecting the truly influential features for high dimensional setting of p>n, it may neglect some important features that are jointly associate with the response but have weak individual effects. Furthermore, as a marginal approach, it may rank highly some unimportant SNPs that are spuriously correlated with the response due to their strong collinearity with other influential features [26, 37]. To overcome these shortcomings, we use the iterative process to address possible complex situations of SNPs that can exist.

The main difference between TC-SIS and ITC-SIS is that TC-SIS finalizes the first d members of by only one step while IDC builds up gradually with several steps [37], i.e. , with d=d1+d2+…+d, where stands for the SNPs selected at i step and d is the number of SNPs for each set , for i=1,…,k. The main idea of ITC-SIS is to iteratively adjust residuals obtained from regressing all remaining SNPs onto the selected ones contained in . Regressing unselected on selected, and adjusting residuals, effectively breaks down original complex correlation structure among SNPs. The iterative steps of ITC-SIS can be summarized as:

Step 1: Use to rank all SNPs based on their individual trend correlations with the response, and then input the first d1 members into , where d1

Step 2: Define , where is the complement set of . Then choose the second d2 members into using TC-SIS to rank all candidates of X for Y, where d1+d2≤d.

Step 3: repeat step 2 until the size of reaches the pre-specified number d.

See more details of the iterative process from Zhong et al. [37].

We refer the reader to an example of a single gene with strong SNP.

Results

In this section, we assess the performance of TC-SIS by four empirical Monte Carlo simulation studies under various designs, and also a real data analysis examining the genome-wide association studies on PCOS affection status. We evaluate the performance of the screening procedures through the following three criteria [24]:

Average Minimum Model Sizes : The average of the minimum number of features that are required by each screening procedure to select all truly influential features across all simulation replicates. The closer to the true model size for the estimated is, the better the screening procedure is determined to be.

Individual Success Rates : The proportion of each truly influential feature is correctly selected by the screening method within the threshold d across all simulation replicates. This requires that the screening score of each truly influential feature ranks within the top d among all p features.

Simultaneous Success Rates (P): The proportion of replicates in which all of the truly influential features are simultaneously selected by the screening method within the threshold d. This requires that the screening scores of all truly influential features rank within the top d among all p features. The closer to one that this proportion is, the better the screening procedure is determined to be.

3.1 Simulation study 1

In this Simulation Study we directly adopt a published real genome-wide association data collected from rice accessions, including 36,901 SNPs and 272 samples [38]. We set the first five SNPs as the truly influential ones, and notice that the correlations among them are complex, ranging from 0.19 (slightly correlated) to 1 (perfectly correlated). The remaining 36,896 SNPs serve as confounding noise, representing a complicated genome simulation setting. We design a joint effect of these five loci by generating the response from a multiple logistic regression model as where the residual term ε is generated from N(0,1). The coefficients of the five influential SNPs, β(i=1,...5), are randomly selected from a mixed Gaussian distribution:

In Simulation Study 1, we replicate each simulation 100 times, and compare four methods: DC-SIS, MMLE-SIS, TC-SIS, and ITC-SIS. For ITC-SIS, d1 is set to be 6, k=2, and d2=d−d1. The results are summarized in Table 1. Given each of the same thresholds, TC-SIS achieves higher individual success rates than DC-SIS and MMLE-SIS. However, the simultaneous success rates of all the first three methods are all very low (close to zero) because they are trapped by a couple of influential loci each having very weak individual effect but associating with the response by joint effects with other loci. Therefore, a large amount of confounding SNPs that are not actually associated with the response but appear to be important because of their high LDs with the other loci act as a role to confuse the individual/marginal approaches (the first three) to include 13000 SNPs on average to locate all of the five true loci. Compared to these three individual approaches, ITC-SIS requires an average model size of only 31.69 to simultaneously select all of the five truly influential loci from 36,901 SNP candidates across the 100 simulation replicates. It is a striking improvement because 31.69 is only 2 thousandth (0.002) of the model size needed by the first three approaches. It indicate that the iterative process is very effective in successfully detecting the true multi-loci without being trapped by spurious associations caused by LD. The running time of one replicate for ITC-SIS is around 11 seconds on a MacBook Pro with 2.2 GHz Intel Core i7 and 16GB RAM.

Table 1

Success rates of four feature screening approaches in selecting each and all truly influential feature X within thresholds d=20,40,60 for Simulation Study 1

d=20
	Pa	|ℳ^|	P_X_1	P_X_2	P_X_3	P_X_4	P_X_5	Run Time
DC-SIS	0	20927.89	0.45	0.46	0.27	0.60	0.48
MMLE-SIS	0	13982.14	0.77	0.78	0.25	0.39	0.78
TC-SIS	0.03	28202.48	0.99	0.67	0.29	0.77	0.45	5.042s
ITC-SIS	0.22	31.69	1.00	0.99	0.47	0.76	0.99	11.025s
d=40
	Pa	|ℳ^|	P_X_1	P_X_2	P_X_3	P_X_4	P_X_5	Run Time
DC-SIS	0	20927.89	0.54	0.55	0.30	0.61	0.56
MMLE-SIS	0	13982.14	0.79	0.79	0.27	0.39	0.79
TC-SIS	0.03	28202.48	0.99	0.67	0.31	0.77	0.45	5.042s
ITC-SIS	0.89	31.69	1.00	1.00	0.95	0.94	1.00	11.025s
d=60
	Pa	|ℳ^|	P_X_1	P_X_2	P_X_3	P_X_4	P_X_5	Run Time
DC-SIS	0	20927.89	0.61	0.62	0.31	0.62	0.63
MMLE-SIS	0	13982.14	0.79	0.79	0.27	0.39	0.79
TC-SIS	0.03	28202.48	0.99	0.67	0.31	0.78	0.45	5.042s
ITC-SIS	0.94	31.69	1.00	1.00	0.97	0.97	1.00	11.025s

For each of the following three simulation designs, we fix the sample size, n, to be 200 and set the number of features, p, to be 5,000. We replicate each simulation 500 times and compared four approaches, MMLE-SIS, DC-SIS, PC-SIS, and TC-SIS. To be fair, we compared only these four marginal approaches without applying iterative process for any of them.

3.2 Simulation study 2

Each observation of the response, Y, will be generated by a Bernoulli process with where p∼Unif(0.05,0.95) is chosen anew for each replicate of the simulation. We design the first ten features to be truly associated with the response Y, i.e., S={1,…,10}. Similar to Example 1 of [2], we generate these first ten features as with the values of π being given by Table 2. This means that each causative X will take on values of 0, 1, or 2 (representative of three ordinal levels, with 0≺1≺2). For any j∉S, we generate X∼Binomial(2,π) with π∼Unif(0.05,0.95). The value of π is chosen anew with each replicate of the simulation. This means that these non-causative features will have no association with Y.

Table 2

Values of π used to simulate data in Simulation Study 2

	πm1	πm2	πm3	πm4	πm5	πm6	πm7	πm8	πm9	πm,10
Y=0	0.3	0.4	0.6	0.7	0.2	0.4	0.3	0.8	0.4	0.2
Y=1	0.6	0.1	0.1	0.4	0.8	0.7	0.9	0.2	0.7	0.6

The results are summarized in Table 3. For this simulation, TC-SIS results in the smallest average model size of 54.674 to contain all the ten truly influential features, which is ten features less than the next closest method (DC-SIS for 64.990). The average minimum model sizes of PC-SIS and MMLE-SIS nearly double or triple that required by TC-SIS, respectively. In the case of TC-SIS versus MMLE-SIS, the individual success rates are at times nearly fourfold more favorable towards our method.

Table 3

Success rates of four feature screening approaches in selecting each truly influential feature X within thresholds d=15, respectively for Simulation Study 2

d=15
	|ℳ^|	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10
MMLE-SIS	150.340	0.384	0.746	0.756	0.404	0.822	0.844	0.354	0.742	0.320	0.400
DC-SIS	64.990	0.900	0.984	0.990	0.898	0.998	0.998	0.894	0.984	0.888	0.894
PC-SIS	93.018	0.862	0.974	0.980	0.864	0.994	0.998	0.860	0.966	0.854	0.862
TC-SIS	54.674	0.916	0.988	0.994	0.922	1.000	0.998	0.912	0.982	0.904	0.908

3.3 Simulation study 3

Inspired by the concept of discretization of a continuous random variable as found in Example 3 of [2], we connect the influential features with the response via an indirect way. Similar to Simulation Study 2, we generate the response Y from a Bernoulli process with where p∼unif(0.05,0.95) is again chosen anew for each replicate of the simulation. Given Y=m, we generate a latent variable Z independently distributed as N(Y,1) for the first ten truly influential features j∈S. The first ten influential features X are then discretized from Z based on the cutoffs (κ,κ) listed in Table 4 as:

Table 4

Values of (κ,κ) used to simulate data in Simulation Study 3

	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10
κLj	0	0	0.2	0	-0.2	0.2	0	0.1	-0.2	0.2
κUj	0.7	1	0.8	0.9	1.2	1	1	1	1.2	0.8

These cutoffs in Table 4 are set to establish weaker associations between the response Y and each of the influential feature to increase the difficulty level in recognizing the true model. It should be noted that this method of generating the truly influential features results in a trend association between each of the truly influential feature and response: namely, lower values of X are associated with Y=0 and higher values of X are associated with Y=1. For any j∉S, we generate X∼Binomial(2,π) with π∼Unif(0.05,0.95).

The results are summarized in Table 5. TC-SIS here results in the smallest average minimum model size of 112.627 to get the true model. This leads us to the conclusion that TC-SIS does a better job than other approaches at avoiding ballooning models. Of especial note here, MMLE-SIS fails on average to produce a selected model smaller than the sample size of n=200 and its success rates are at many times less than 0.1 (versus 0.8 of TC-SIS). These results demonstrate the capability of TC-SIS to obtain excellent results when trend correlation exists between the feature and response. The performance of DC-SIS is comparable in the success rates but at the cost of a relatively larger model.

Table 5

Success rates of four feature screening approaches in selecting each truly influential feature X within thresholds d=15, respectively for Simulation Study 3

d=15
	|ℳ^|	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10
MMLE-SIS	508.672	0.072	0.060	0.066	0.078	0.554	0.054	0.388	0.098	0.032	0.204
DC-SIS	125.258	0.876	0.884	0.886	0.904	0.886	0.880	0.880	0.910	0.878	0.906
PC-SIS	171.829	0.820	0.806	0.816	0.842	0.876	0.806	0.876	0.830	0.820	0.858
TC-SIS	112.627	0.876	0.876	0.882	0.904	0.924	0.874	0.920	0.906	0.878	0.900

3.4 Simulation study 4

In this simulation, we generate the sample data using a logistic regression model, which is the fundamental basis of MMLE-SIS. We first generate each feature X (1≤j≤p) by uniformly sampling from the set {0,1,2} with equal probability. We then connect a binary response Y with the first five features by letting and using

to sample the binary response. The coefficients are as given in Table 6. The results are summarized in Table 7. For this example, we once again obtain a smaller required average minimum model size than DC-SIS and PC-SIS (46.470 for DC-SIS, 93.270 for PC-SIS, as compared to 41.976 for TC-SIS). Unlike the first two examples, MMLE-SIS recoups its earlier collapses and matches TC-SIS nearly perfectly in both success rates and average minimum model sizes in Example 3. Although it is the specific strength of MMLE-SIS to handle logistic regression model, TC-SIS still produces results abreast with that of MMLE-SIS. In addition, it should be noted that since MMLE-SIS requires solving an optimization problem to produce its screening statistics, TC-SIS is significantly faster in computational run time. Thus, when run time is an issue, we suggest the use of TC-SIS over MMLE-SIS, even when the logistic regression model holds.

Table 6

Values of used to simulate data in Simulation Study 4

	β_X_1	β_X_2	β_X_3	β_X_4	β_X_5
Xj=0	0	-5	2	-6	1
Xj=1	3	-3	4	-4	3
Xj=2	5	-1	6	-2	5

Table 7

Success rates of four feature screening approaches in selecting each truly influential feature X within thresholds d=15, respectively for Simulation Study 4

d=15
	|ℳ^|	X1	X2	X3	X4	X5
MMLE-SIS	41.934	1.000	0.856	0.868	0.842	0.870
DC-SIS	46.470	1.000	0.860	0.850	0.838	0.866
PC-SIS	93.270	1.000	0.794	0.758	0.778	0.790
TC-SIS	41.976	1.000	0.860	0.858	0.842	0.862

3.5 Real data analyses

We apply the proposed ITC-SIS screening procedure to a clinical dataset pertaining to the genome-wide association studies on the PCOS affection status (dbGaP Study Accession: https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000368.v1.p1). This data consists of 4,099 subjects (1,043 cases and 3,056 controls) and 731,442 SNPs. The goal of this analysis is to identify the most influential susceptibility loci that affect PCOS status for European Caucasian population. The response for this data is PCOS affection status (binary) and the features are the encoded SNP genotype values (categorical), which exactly represents a problem that the ITC-SIS is originally motivated.

We removed rare alleles that have minor allele frequency (MAF) less than 0.1 before performing the ITC-SIS. To determine the optimal value d1 for the fist iterative process, we investigate each value among a set of d1=1,2,…,3[n4/5/log(n4/5)]=351 and check the mean square prediction error (MSPE). The prediction is assessed by a cross validation process, with 75% of the observed data for training and 25% for testing. As shown in Fig. 1, as the model size increases, the MSPE first drops and then stays flat. We choose d1 as the minimum model size whose MSPE falls into one standard deviation plus the minimum MSPE. As demonstrated in Fig. 1, we select d1=175 SNPs in the first iteration, d2=351−175=176 SNPs in the second iteration, and number of iterations k=2. The remaining 731,091 SNPs are filtered out as noise.

Fig. 1

The MSPE used to select d1 in Real Data Analyses

After the completion of ITC-SIS, high dimensionality is not an issue any more (n=4099,d=351). Then we apply the multiple logistic regression model to estimate the multi-loci joint effects of these promising candidates, which is able to assess the significance level of each SNP through p-values. We compare the results of multiple logistic regression after integrating it with each of the DC-SIS, TC-SIS, and ITC-SIS screening process (see Table 8). Three model selection criteria are used: model size, Akaike’s Information Criterion (AIC), and misclassification rate, which is computed as the percentage of incorrectly predicted affection status after applying the multiple logistic regression model to fit the selected 351 SNPs that are selected by each screening procedure.

Table 8

Model Selection of Three approaches Applied in Real Data Analyses

Two-stage Method	Model size	AIC	Misclassification Rate
DC-SIS + Multiple Logistic Regression	70	4188.33	21.74%
TC-SIS + Multiple Logistic Regression	86	3601.02	19.51%
ITC-SIS (d1=175) + Multiple Logistic Regression	88	3581.96	19.27%

As expected, ITC-SIS + multiple logistic regression yields the best model with the smallest misclassification rate and the smallest AIC. It suggests 88 influential SNPs that was highlighted as red triangle in Fig. 2. Figure 2 is very different from the traditional Manhattan plot that was obtained by single-SNP approaches from several aspects: 1) It is much less dense because the iterative feature screening process dramatically shrinks noise SNPs into zero that makes 99.9% of SNPs disappeared from current plot. The noise SNPs built up a very tall and dense base in traditional Manhattan plot. 2) Unlike traditional Manhattan plot, the vertical axis demonstrates ITC-SIS scores instead of p-values. As a two-stage approach that explore the complex structure of the data, individual p-value will not be meaningful in this plot. 3) It separates important SNPs from noise SNPs in a much striking way. Specifically, noise SNPs are in the bottom line that makes the selected SNPs substantially stand out. 4) It revolutionizes the traditional selection rule that only select the significance from the top of the traditional Manhattan plot. You may wonder why several SNPs with very small scores are selected but others with higher scores are not selected. It actually reflects the joint effects of multi-loci and confounding issue of LD described in early sections of this article. Specifically, PCOS is a complex disease that are affected by multi-loci each having small individual effects, meanwhile many noise SNPs may have strong signals because they are in strong multicollinearity with other truly influential SNPs but they actually do not directly associated with the disease. To show the data-driven nature of our proposed approach, we refer the readers to a publication that also used iterative feature screening for GWAS data but worked on a single-gene trait with strong individual effects. As you can see, their SNPs were selected based on scores locating on the top of the Manhattan plot [15].

Fig. 2

The Manhattan plot of the ITC scores in Real Data Analyses

There are over 50 genes being located from these 88 informative SNPs. Additional file 1: Table 1 in the Supplement file summarizes the estimated ITC score, coefficient, p-value, corresponding gene name, Allele type, and detailed position for each of the 53 selected influential SNPs that could locate nearby genes. In addition to confirming many genes (FSHR, LHCGR, C9orf3, RPS26, RAB5B, SUOX, ERBB3, TOX3, ApoB, ROBO2, NEIL2) that were reported to be directly associated with PCOS, we also detect several new genes. Specifically, we find that the SNP rs7559066 located at Chr 2 lies within the FSHR gene and the SNP rs12235316 located at Chr 9 lies within the C9orf3 gene. The FSHR and C9orf3 gene have been reported by numerous studies to be strongly associated with PCOS in women and erectile dysfunction in men for both Han Chinese population [39] and European Caucasian population [40], which appear in individuals who have either inadequate or excessive amounts of sexual hormones. Adolescent girls with obesity and PCOS were found to have elevated fasting and postprandial plasma TG and ApoB-lipoprotein remnants [41]. ROBO2 gene were found differentially expressed between obese women without PCOS and obese women with PCOS [42]. NEIL2 gene may help identify pathways that link specific PCOS related traits with greater metabolic risk [43]. ERBB3 is T2D candidate gene, implicated in the process of female gamete generation and determining function of antigen-presenting cells [39].

The new genes detected from this data analyses could be found in Table 1 of the supplement file. To name a few, TACR1 and GASK1A have broad expression in endometrium. LTBP2 has broad expression in ovary (RPKM 21.5). NR2C2 encodes a protein that belongs to the nuclear hormone receptor family [43]. MSH6 and BRCC3 were found to be relevant for PCOS related phenotypes by a new protein-protein interaction network analysis [44]. NR2C2 encodes a protein that belongs to the nuclear hormone receptor family functioning in development, cellular differentiation and homeostasis. In summary this analysis confirm many genes that were reported to be associated with PCOS and also locate several new genes that are related to endometrium, hormones, organ growth, and cell division. Their functions in PCOS need further investigations by molecular and functional genetics.

Discussion

In this paper, we propose a new feature screening procedure using trend correlation, whose finite performance is demonstrated via performing multiple simulation studies with various designs and also comparing with three other relevant extant approaches. We furthermore illustrate the performance of TC-SIS through real data analyses pertaining to genome-wide association studies on PCOS disease. We establish the strong screening consistency for this procedure when the number of features diverges exponentially with respect to the sample size. Strong screening consistency is much harder to achieve than the sure screening property, as it guarantees that not only the selected model contains the true model, but also that the selected model equals to the true model asymptotically.

The proposed TC-SIS method can be easily extended to a categorical response having greater than two levels if needed; however, we only consider binary Y here because this allows for some simplification of our notation and proofs. It has been noted that the choice of a threshold d, the number of SNPs to keep, is of importance in feature screening literature. Several methods have been proposed to determine such a threshold, e.g. [2, 29, 45, 46]. We follow the rule of thumb proposed by Liu et al. [47], and set the cutoff of model size as multiplier of d=[n4/5/log(n4/5)]. In the simulation study 1 (a harder case), we test the cutoffs d=20, 2d=40, and 3d=60 for n=272. We set d=15 when n=200 for all other simulation studies. In the real data example, we choose a model size 3d=351 when n=4099 to avoid missing influential candidate from the beginning, and then test the significance of these candidate SNPs by a joint model that performs well if high dimensionality is not an issue.

In addition to the general association detected by other methods, TC-SIS excels in exploring the trend association between the response and the features, e.g., larger feature values tends to be associated with larger (or conversely, smaller) response values in certain practices. Another appealing advantage that we observe from the simulation studies is the relative stability of TC-SIS (compare to other methods) in the face of a potentially large unbalance in the number of positive (Y=1) responses. TC-SIS assumes milder conditions than other approaches in that it neither requires any regression model structure nor assumes any specific distribution of the data.

Conclusion

Detecting important multi-loci that are associated with the complex disease is challenging because each locus may have weak effect. The ITC-SIS following by a multiple regression model serves as a good alternative method to detect disease susceptibility loci for clinic genomic data. It confirms around ten genes that were reported to be associated with PCOS and also detects many new genes after scanning a high dimensional set of SNPs.

Additional file 1 In the Supplementary Materials we present in full the proofs for Theorems 1 and 2 given at “Theoretical properties” section of the main text, and additional results from real data analyses.