Simultaneous selection of multiple important single nucleotide polymorphisms in familial genome wide association studies data

We propose a resampling-based fast variable selection technique for detecting relevant single nucleotide polymorphisms (SNP) in a multi-marker mixed effect model. Due to computational complexity, current practice primarily involves testing the effect of one SNP at a time, commonly termed as ‘single SNP association analysis’. Joint modeling of genetic variants within a gene or pathway may have better power to detect associated genetic variants, especially the ones with weak effects. In this paper, we propose a computationally efficient model selection approach—based on the e-values framework—for single SNP detection in families while utilizing information on multiple SNPs simultaneously. To overcome computational bottleneck of traditional model selection methods, our method trains one single model, and utilizes a fast and scalable bootstrap procedure. We illustrate through numerical studies that our proposed method is more effective in detecting SNPs associated with a trait than either single-marker analysis using family data or model selection methods that ignore the familial dependency structure. Further, we perform gene-level analysis in Minnesota Center for Twin and Family Research (MCTFR) dataset using our method to detect several SNPs using this that have been implicated to be associated with alcohol consumption.


Introduction
Genome Wide Association Studies (GWAS) have identified a large number of genetic variants associated with complex diseases [4,50].The advent of economical high-throughput genotyping technology enables researchers to scan the genome with millions of Single Nucleotide Polymorphism (SNP)-s, and improvements in computational efficiency in analysis techniques has facilitated parsing through this huge amount of data to detect significant associations [44].However, detecting small effects of individual SNPs requires large sample size [32].For quantitative behavioral traits such as alcohol consumption, drug abuse, anorexia and depression, variation in genetic effects due to environmental heterogeneity brings in additional noise, further amplifying the issue.This is one of the motivations of performing GWAS on families instead of unrelated individuals, through which the environmental variation can be reduced [2,35,45].However, association analysis of multiple SNPs while using dependent data with a familial structure and large sample sizes can be computationally very challenging.Thus single SNP association analysis is the standard tool for detecting SNPs, and most family studies tend to have smaller sample size.The MCTFR Study [35] with genome-wide data on identical twins, non-identical twins, biological offspring, adoptees serve as the motivation for our methodology development in this paper.
A downside of family-based single-SNP methods-such as GRAMMAR [1] and the association test of Chen and Abecasis [7]-is that they do not take into account shared environment effects within families.They assume that phenotypic similarity among individuals in a family is entirely due to their genetic similarity and not due to the effect of shared environment.As a result, they tend to lose power when analyzing data where shared environmental effects explain a substantial proportion of the total phenotypic variation (see [11,34] for examples).The RFGLS method proposed by Li et al. [28] does take into account genetic and environmental sources of familial similarity and provides fast inference through a rapid approximation of SNP-specific coefficients from a mixed effect model.However it is only able to 1 arXiv:1802.01141v3[stat.AP] 21 May 2023 handle single SNPs at a time.
Single-SNP methods are less effective in detecting SNPs with weak signals [32].This is limiting in situations where multiple SNPs are jointly associated with the phenotype [22,40,52].Several methods of multi-SNP analysis have been proposed as alternatives.The kernel based association tests [6,19,39,40] are prominent among such techniques.However, all such methods test for whether a group of SNPs is associated with the phenotype of interest as a whole, and do not prioritize within that group to detect the individual SNPs primarily associated with the trait.One way to solve this problem is to perform model selection.The methods of Frommelet et al. [14] and Zhang et al. [53] take this approach, and perform SNP selection from a multi-SNP model on GWAS data from unrelated individuals.However, they rely on fitting models corresponding to multiple predictor sets, hence are computationally very intensive to implement in a linear mixed effect framework for modeling familial data.
In this paper we propose a fast and scalable model selection technique that fits a single model to a familial dataset, and aims to identify genetic variants with weak signals that are associated with the outcome through joint modelling of multiple variants.We consider only main effects of the variants, but this can be extended to include higher-order interactions.We achieve this by extending the framework of e-values [31], which we discuss in Section 2.4.There we present the definition of an e-value, discuss some of its properties, and describe how we generalize the e-values for our scenario.Broadly, for any estimation method that provides consistent estimates (at a certain rate relative to the sample size) of the vector of parameters, e-values quantify the proximity of the sampling distribution for a restricted parameter estimate to that of the full model estimate in a regression-like setup.A variable selection algorithm using the e-values has the three simple and generic steps.First, fit the full model, i.e.where all predictor effects are being estimated from the data, and use resampling to estimate its e-value.Second, set an element of the full model coefficient estimate to 0 and get an e-value for that predictor using resampling distribution of previously estimated parameters-repeat this for all predictors.Then finally, select predictors that have e-values below a pre-determined threshold.
The above algorithm offers multiple important benefits in the SNP selection scenario.Unlike other model selection methods, only the full model needs to be computed here.It thus offers the user more flexibility in utilizing a suitable method of estimation for the full model.Our method allows for fitting multi-SNP models, thereby accommodating cases of modelling multiple correlated SNPs or closely located multiple causal SNPs simultaneously.Finally, we use the Generalized Bootstrap (GBS) [5] as our chosen resampling technique.Instead of fitting a separate model for each bootstrap sample, GBS computes bootstrap estimates using Monte-Carlo samples from the resampling distribution as weights, and reusing model objects obtained from the full model.Consequently, the resampling step becomes very fast and parallelizable.
In past literature, VanderWeele and Ding [43] and Vovk and Wang [47] used the term 'e-value' in the contexts of sensitivity analysis and multiple testing, respectively.In comparison, the e-values we use [31] evaluate the relevance of a variable with reference to a statistical model.Going beyond the existing proposal of e-values tied to specific objectives and models, as well as the well-known p-values used for hypothesis testing, this e-value is assumption-lean, covers more generic statistical problemssuch as including dependent data models-and is expandable to numerous applications, including group feature selection, hypothesis testing, and multiple testing.

The MCTFR data
The familial GWAS dataset collected and studied by Minnesota Center for Twin and Family Research (MCTFR) [28,34,35] consists of samples from three longitudinal studies conducted by the MCTFR: (1) the Minnesota Twin Family Study (MTFS) [18] that covers twins and their parents, (2) the Sibling Interaction and Behavior Study (SIBS) [33] that includes adopted and biological sibling pairs and their 2/29 rearing parents, and (3) the enrichment study [23] that extended the MTFS by oversampling 11 year old twins who are highly likely to develop substance abuse.While 9827 individuals completed the phenotypic assessments for participation in the study, after several steps of screening [35] the genotype datafrom 7605 Caucasian individuals clustered in 2151 nuclear families were included in our analysis.This consisted of 1109 families where the children are identical twins, 577 families with non-identical twins, 210 families with two adopted children, 162 families with two non-twin siblings, and 93 families where one child is adopted while the other is the biological child of the parents.
DNA samples collected from the subjects were analyzed using Illumina's Human660W-Quad Array, and after standard quality control steps [35], 527,829 SNPs were retained.Covariates for each sample included age, sex, birth year, generation (parent or offspring), as well as the two-way interactions generation x age, generation x sex, and generation x birth year.Five quantitative phenotypes measuring substance use disorders were studied in this GWAS: (1) Nicotine dependence, (2) Alcohol consumption, (3) Alcohol dependence, (4) Illegal drug usage, and (5) Behavioral disinhibition.The response variables corresponding to these phenotypes are derived from questionnaires using a hierarchical approach based on factor analysis [16].
A detailed description of the data is available in Miller et al. [35].Several studies reported SNPs associated with phenotypes collected in MCTFR study [8,28,34].Li et al. [28] used RFGLS to detect association between height and genetic variants through single-SNP analysis, while McGue et al. [34] used the same method to study SNPs influencing the development of all five indicators of behavioral disinhibition mentioned above.Irons [20] focused on the effect of several factors affecting alcohol use in the study population, namely the effects of polymorphisms in the ALDH2 gene and the GABA system genes, as well as the effect of early exposure to alcohols as adolescents to adult outcomes.Finally Coombes et al. [8] used a bootstrap-based combination test and a sequential score test to evaluate gene-environment interactions for alcohol consumption.

Consents and Approvals
Data were collected through the Minnesota Center for Twin and Family Research (MCTFR).All University of Minnesota and National Institute of Health (NIH) guidelines for human subjects research were followed in the collection and processing of the data.The protocol was approved by the Institutional Review Board (IRB) at the University of Minnesota (protocol # 0303M45703).Participants aged 18 years and older completed informed consent, while consent was obtained from at least one parent for those participants younger than 18 and the minor participant also assented to participate.

Statistical model
We use a Linear Mixed Model (LMM) with three variance components accounting for several potential sources of variation to model effect of SNPs behind a quantitative phenotype.This is known as ACE model in the literature [24].While the-state-of-the-art focuses on detection of a single variant at a time, we will incorporate all SNPs genotyped within a gene (or group of genes in some cases) as set of fixed effects in a single model.
Following standard protocol for family-based GWAS [7,28,34], we assume a data setting of nuclear pedigrees, i.e. that the data consists of observations from individuals of multiple genetically unrelated families, with individuals within a family potentially sharing genetic material.Suppose there are m such families in total, with the i th pedigree containing n i individuals, and the total number of individuals is n = m i=1 n i .Denote by y i = (y i1 , . . ., y in i ) T the quantitative trait values for individuals in that pedigree, while the matrix G i ∈ R n i ×pg contains their genotypes for a number of SNPs.Let C i ∈ R n i ×p denote the data on p covariates for individuals in the pedigree i.Given these, we consider the following model.
with α the intercept term, β g and β c fixed coefficient terms corresponding to the multiple SNPs and covariates, respectively, and i ∼ N n i (0, V i ) the random error term.To account for the within-family dependency structure, we break up the random error variance into three independent components: The three components of V i in (2.2) model different sources of random variations that can affect the quantitative trait values for individuals in the i th pedigree.The first component above is a within-family random effect term to account for shared polygenic effects.The proportion of of genetic material shared between pairs of individuals in a family is represented by elements of the matrix Φ i .Its (s, t) th element represents two times the kinship coefficient, which is the probability that two alleles, one randomly chosen from individual s in pedigree i and the other from individual t, are 'identical by descent', i.e. come from same common ancestor [24].Following basic probability, the kinship coefficient of a parent-child pair is 1/4, a full sibling pair or non-identical (or dizygous = DZ) twins is 1/4, and for identical (or monozygous = MZ) twins is 1/2 in a nuclear pedigree.Following this, we can construct the Φ i matrices for different types of families: for families with parents (indices 1 and 2) and MZ twins, DZ twins, or two adopted children (indices 3 and 4), respectively.
The second variance component σ 2 c 11 T in (2.2) accounts for shared environmental effect within each pedigree.Traits of each individual in the pedigree are affected by the same amount-a single random draw from N (0, σ 2 c )-of random variation.The third term in (2.2) quantifies other sources of variation unique to each individual.

Feature Selection with e-values
We extend the recently-proposed framework of e-values [31] to select important SNPs in the above genelevel, multi-SNP statistical model.In a general modelling situation where one needs to estimate a set of parameters θ ∈ R d from data with sample size n, a statistical model corresponds to a subset of the full parameter space.In other words, the estimable index set of θ, say S ⊆ {1, . . ., d} specifies a model.The other indices are set at constant values-typically in model selection literature the constants are set at 0. Note that we are attempting to select important SNPs as described in Section 2.3, thus in our setting θ ≡ β g , d ≡ p g .
Following the recipe in Majumdar and Chatterjee [31], we obtain coefficient estimates corresponding to model S by simply replacing elements of the 'full model' estimate θ-i.e. the the coefficient estimate with all possible parameters included-at indices not in S: Sampling distribution is defined as the distribution of a parameter estimate, based on the random data samples used to calculate this estimate.We compare sampling distributions of the above model with the full model, i.e. [ θS ] with [ θ] (denoting the distribution of a random variable by [•]).For this comparison, we define an evaluation map function E : R d × Rd → [0, ∞) that measures the relative position of θS with respect to [ θ].Here Rd is the set of probability measures on R d .For any x ∈ R d and [X] ∈ Rd with a positive definite covariance matrix VX, we consider the following evaluations functions in this paper: 3) Here diag(VX) denotes the vector composed of the diagonal entries in VX, and represents elementwise product, so that (diag(VX)) 1/2 is the vector of coordinate-wise standard deviation, and (diag(VX)) −1/2 (x − EX) is a normalized version of x.Data depths [42,55] also constitute a broad class of functions that can be used as evaluation maps-as done by Majumdar and Chatterjee [31].In general, any continuous function that is location and scale invariant, and has a few basic convergence properties is a good choice for the evaluation map function (see conditions E1-E4 in Appendix A.1).

Formulation
Note that the evaluation map function is defined conditional on a fixed value of θS .Since θS itself has a distribution, so does the evaluation map.We define the e-value as any functional of the evaluation map distribution E S ≡ [E( θS , [ θ])] that can act as a measure of comparison between the sampling distributions of θS and θ.For example, Majumdar and Chatterjee [31] took the mean functional of E S (say µ(E S )) as e-value, and showed that it can be used as a model selection criterion.To this end, non-zero indices (say S 0 ) of the true parameter vector θ 0 can be recovered through a fast algorithm that has these generic steps: Algorithm 1: Model Selection with e-values 1. Obtain the e-value of the full model, say µ(E * ). 2. For the j th predictor, j = 1, . . ., p, consider the model with the j th coefficient of θ replaced by 0, and obtain its e-value-denoted by µ(E −j ).
3. Collect the predictors for which µ(E −j ) < µ(E * ): denote it by Ŝ0 .This is the estimated set of non-zero coefficients in θ.
As n → ∞, the above algorithm provides consistent model selection, i.e.P( Ŝ0 = S 0 ) → 1.In practice we only have one dataset, so it is not possible to access the true sampling distribution of θ and θS to do the above.To this end, we use a fast bootstrap algorithm, called Generalized Bootstrap (GBS) [5], to obtain approximations of the sampling distributions [ θS ], [ θ], the evaluation map distributions, and the e-values.GBS is dependent of a tuning parameter τ n that represents the standard deviation of the synthetic noise introduced by the bootstrap procedure.Intermediate values of τ n , such that τ n /n → ∞, result in model selection consistency as described above.

Quantile e-values
When true signals are weak, the above method of variable selection leads to very conservative estimates of non-zero coefficient indices, i.e. a large number of false positives in a sample setting.This happens because the true values leave-one-covariate-out e-values for variables that correspond to small but nonzero coefficients in θ 0 (hence weak signal) fall too close to the full model e-value.Consequently, when these e-values are estimated from randomly sampled data, simply by random chance their values can be slightly less than the full model e-value estimate.Figure 1 demonstrates this phenomenon in our setup, where we would like to estimate non-zero elements of the fixed effect coefficient vector β g in the model (2.1), i.e. β g ≡ θ.
Here we analyze data on 250 families with monozygotic twins, each individual being genotyped for 50 SNPs.Four of these 50 SNPs are causal: each having a heritability of h/6% with respect to the total error variation present.The four panels show density plots of E −j for j = 1, . . ., p, as well as E * : estimated based on resampling schemes with four different values of the standard deviation parameter s ≡ s n = τ n / √ n.Focusing on where the central regions of the evaluation map distributions are, we notice that for smaller values of s there is quite a bit of overlap along the bootstrap estimates of E −j for causal vs. non-causal SNPs.On the other hand, for large values of s all the density plots become essentially the same as the full model.
However, notice that the evaluation map distributions for non-zero vs. zero indices have different tail behaviors at smaller values of s.In the Appendix we show that the means and tail quantiles for E −j and E * asymptotically converge to different limits (Theorems A.1 and A.2), and these limits are well-approximated with a GBS scheme having small standard deviation s (Theorem A.3).A potential reason for the different tail behaviors we see above is that the convergence at tail quantiles happens at a faster rate than convergence at the means.Consequently, instead of comparing means of the distributions, comparing a suitable tail quantile across the distributions is more likely to provide a better separation of non-zero vs. zero indices.For this reason we use tail quantiles as e-values.
When the q th quantile (denoted by c q , q ∈ (0, 1)) is taken as the e-value instead of the mean, we set a lower detection threshold than the same functional on the full model, i.e. choose all j such that to be included in the model.The optimal choice of q and t depends on factors such as specifications of the statistical model, sample size, and degree of sparsity of parameters in the data generating process.We demonstrate this point through our experiments in Section 3.For q, we take the conservative route by only flagging a SNP as 'detected' if Eq. (2.5) holds for all q ∈ {0.5, 0.6, 0.7, 0.8, 0.9}.This approach leads to a tradeoff between the true positive and true negative SNP detections rates for different values of t.We demonstrate this fact through synthetic data experiments (Section 3.1), and choose the best t that minimizes prediction error on a holdout sample in the MCTFR data analysis (Section B).

Experiments
We now evaluate the performance of the above formulation of quantile e-values in through on synthetic data, as well as the MCTFR Twin Studies dataset. 6/29

Synthetic Data
Consider the model in (2.1) with no environmental covariates and familes with MZ twins.We take a total of p g = 50 SNPs, and generate the SNP matrices G i in correlated blocks of 6, 4 ,6, 4 and 30 to simulate correlation among SNPs in the genome.We set the correlation between two SNPs inside a block at 0.7, and consider the blocks to be uncorrelated.For each parent we generate two independent vectors of length 50 with the above correlation structure, and entries within each block being 0 or 1 following Bernoulli distributions with probabilities 0.2, 0.4, 0.4, 0.25 and 0.25 (Minor Allele Frequency or MAF) for SNPs in the 5 blocks, respectively.The genotype of a person is then determined by taking the sum of these two vectors: thus entries in G i can take the values 0, 1 or 2. Finally we set the common genotype of the twins by randomly choosing one allele vector from each of the parents and taking their sum.We repeat the above process for m = 250 families.In GWAS generally each associated SNP explains only a small proportion of the overall variability of the trait.To reflect this in our simulation setup, we assume that the first entries in each of the first four blocks above are causal, and each of them explains h/(σ 2 a + σ 2 c + σ 2 e )% of the overall variability.The term h is known as the heritability of the corresponding SNP.The value of the non-zero coefficient in k-th block: k = 1, ..., 4, say β k is calculated using the formula: We fix the following values for the error variance components: σ 2 a = 4, σ 2 c = 1, σ 2 e = 1, and generate pedigreewise response vectors y 1 , . . ., y 250 using the above setup.To consider different SNP effect sizes, we repeat the above setup for h ∈ {10, 7, 5, 3, 2, 1, 0}, generating 1000 datasets for each value of h.

Competing methods
We compare our e-value based approach using the evaluation maps E 1 and E 2 in (A.2) with two groups of methods: (1) Model selection on linear model: Here we ignore the dependency structure within families by training linear models on the simulated data and selecting SNPs with non-zero effects by backward deletion using a modification of the BIC called mBIC2.This has been showed to give better results than single-SNP analysis in a GWAS with unrelated individuals [14] and provides approximate False Discovery Rate (FDR) control [3].
(2) Single-marker mixed model: We train single-SNP versions of (2.1) using a fast approximation of the Generalized Least Squares procedure (named Rapid Feasible Generalized Least Squares or RFGLS [28]), obtain marginal p-values from corresponding t-tests and use two methods to select significant SNPs: the Benjamini-Hochberg (BH) procedure, as well as the Local FDR method [12] (LFDR).
For mBIC2 and BH, we choose a conservative FDR level of 0.05 guided by choices in existing work [25,41].Higher FDR values lead to marginal increase in true positive rate but sharp decreases in true negative rate (see definitions of these metrics below).LFDR tends to be much more conservative than global FDR procedures, so we repeated its simulations for a range of FDR values in {0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5}.Table 1 shows the performance metrics at FDR = 0.4, the highest value where the true negative rate is at least as much as the lowest true negative across all settings of e-values as described below.
With the e-value being the q th quantile of the evaluation map distribution, we set the detection threshold value at the t th multiple of q for some 0 < t < 1.This means all indices j such that the q th quantile of the bootstrap approximation of E −j is less than the tq th quantile of the bootstrap approximation of E * get selected as the set of active predictors.To enforce stricter control on the selected set of SNPs we repeat this for q ∈ {0.5, 0.6, 0.7, 0.8, 0.9}, and take the SNPs that get selected for all values of q as the final set of selected SNPs.Guided by empirical experiments, we chose values of t for E 1 and E 2 that clearly demonstrate the tradeoff of TP/TN (or RTP/RTN) rates for the e-values.
Since the above procedure depends on the bootstrap standard deviation parameter s, we repeat the process for s ∈ {0.3, 0.15, . . ., 0.95, 2}, and take as the final estimated set of SNPs the SNP set Ŝt (s) that and finally (4) Relaxed True Negative (RTN), which is the proportion of SNPs in block 5 undetected.We consider the third and fourth metrics to cover situations in which the causal SNP is not detected itself, but highly correlated SNPs with the causal SNP are.This is common in GWAS [14].We average the above proportions over 1000 replications, and repeat the process for two different ranges of t for E 1 and E 2 .

Results
We present the simulation results in table 1.For all heritability values, applying mBIC2 on linear models performs poorly compared to applying RFGLS and then correcting for multiple testing.This is expected because the linear model ignores within-family error components.
Our method works better than the two competing methods for detecting true signals across different values of h: the average TP rate going down slowly than other methods across the majority of choices for t.All the competing methods (mBIC2, RFGLS+BH, RFGLS+LFDR) have very high true negative detection rates, which is matched by our method for higher values of q.Since all reduced model distributions reside on the left of the full model distribution, we expect the variable selection process to turn more conservative at lower values of t.This effect is more noticeable for lower q, indicating that the right tails of evaluation map distributions are more useful for this purpose.Finally for h = 0, we report only TN and RTN values since no signals should ideally be detected.Note also the fact that we report the performance metrics for E 1 , E 2 considering a very conservative selection process: we only mark a SNP j as 'detected' if c q (E −j ) < tc q (E * ) for all q ∈ {0.5, 0.6, 0.7, 0.8, 0.9}.We experimentally observed that relaxing this condition leads to higher strict TP rates but lower strict TN.
RTP performances for all methods are better than the corresponding TP/TN performances.However, for mBIC2 this seems to be due to detecting SNPs in the first four blocks by chance since for h = 0 its RTN is less than TN.Also E 2 seems to perform slightly better than E 1 , in the sense that it yields a higher TP (or RTP) while having the same TN (or RTN) rates.
Among competing methods, it is interesting to notice that LFDR performs better than the other two for small signal values (h ≤ 3), but worse at higher h.This speaks to the strengths of LFDR in lowsignal situations.On the other hand, LFDR is calculated using density estimates of the null and non-null statistic distributions.Since there are only 50 SNPs in our simulation setting (and even less in the real data setting), the resulting instability is a potential reason for its low performance at high values of h.
For model selection we use E 2 as the evaluation function because of its slighty better performance in the simulations.For each gene-level model, We train the LMM in (2.1) on 75% of randomly selected families, perform our e-values procedure for s = 0.2, 0.4, . . ., 2.8, 3, t = 0.1, 0.15, . . ., 0.75, 0.8; and select the set of SNPs that minimizes fixed effect prediction error on the data from the other 25% of families over this grid of (s, t).Note that we consider a wider range of t than in the simulations.This is because of the fact that instead of demonstrating the tradeoff of true positive and true negative rates, our objective here is to actually choose a set of SNPs.For the competing methods, we set FDR levels at 0.05 for mBIC2 and RFGLS+BH, while choose the level that minimizes fixed effect prediction error on the same holdout data as above for RFGLS+LFDR.As seen in Table 2, our e-value based technique detects a much higher number of SNPs than the two competing methods.Our method selects all but one SNP in the genes ALDH2 and COMT.These are small genes of size 50kb and 30kb, respectively, thus SNPs within them have more chance of being in high Linkage Disequilibrium (LD).On the other hand, it does not select any SNPs in SLC6A4 and DRD2.Variants of these genes are known to interact with each other and are jointly associated with multiple behavioral disorders [21,48].
A number of SNPs we detect (or SNPs situated close to them) have known associations with alcoholrelated behavioral disorders.We summarize this in Table 3. Prominent among them are rs1808851 and rs279856 in the GABRA2 gene, which are at perfect LD with rs279858 in the larger, 7188-individual version of our twin studies dataset [20].This SNP is the marker in GABRA2 that is most frequently associated in the literature with alcohol abuse [10], but was not genotyped in our sample.A single SNP RFGLS analysis of the same twin studies data that used Bonferroni correction on marginal p-values missed the SNPs we detect [20]: highlighting the advantage of our approach.In the Appendix, we give a 10/29 gene-wise discussion of associated SNPs (Appendix B), as well as information on all SNPs (Appendix C).
We plot the 90 th quantile e-value estimates in Figures 2, 3 and 4. We obtained gene locations, as well as the locations of coding regions of genes, i.e. exons, inside 6 of these 9 genes from annotation data extracted from the UCSC Genome Browser database [38].Exon locations were not available for OPRM1, CYP2E1 and DRD2.In general, SNPs tend to get selected in groups with neighboring SNPs, which suggests high LD.Also most of the selected SNPs either overlap or in close proximity to the exons, which underline their functional relevance.

Discussion and conclusion
To expand the above approach to a genome-wide scale, we need to incorporate strategies for dealing with the hierarchical structure of pathways and genes: there are only a few genes associated with a quantitative phenotype, which can be further attributed to a small proportion of SNPs inside each gene.To apply the e-values method here, it is plausible to start with an initial screening step to eliminate evidently non-relevant genes.Methods like the grouped Sure Independent Screening [27] and min-P test [49] can be useful here.Following this, in a multi-gene predictor set, there are several possible strategies to select important genes and important SNPs in them.Firstly, one can use a two-stage e-value based procedure.The first stage is same as the method described in this paper, i.e. selecting important SNPs from each gene using multi-SNP models trained on SNPs in that gene.In the second stage, a model will be trained using the aggregated set of SNPs obtained in the first step, and a group selection procedure will be run on this model using e-values.This means dropping groups of predictors (instead of single predictors) from the full model, checking the reduced model e-values, and selecting a SNP group only if dropping it causes the e-value to go below a certain cutoff.Secondly, one can start by selecting important genes using an aggregation method of SNP-trait associations (e.g.Lamparter et al. [26]) and then run the e-value based SNP selection on the set of SNPs within these genes.Thirdly, one can also take the aggregated set of SNPs obtained from running the e-values procedure on gene-level models, then use a fast screening method (e.g.RFGLS) to select a subset of those SNPs.
We plan to study merits and demerits of these strategies and the computational issues associated with them in detail through synthetic studies as well as in the GWAS data from MCTFR.Finally, the current evaluation map based formulation requires the existence of an asymptotic distribution for the full model estimate.We plan to explore alternative formulation of evaluation maps under weaker conditions to bypass this, thus being able to tackle high-dimensional (n < p) situations.
It is important to remember that in a GWAS setting looking for causal factors of polygenic, quantitative traits is a complex problem.Small effect sizes and high amounts of LD-combined with the influence of environmental covariates-can make finding a set of SNPs behind that trait a noisy process.Typically the random effect error-term is earmarked to account for and quantify such heterogeneities at family-level, but how accurate this quantification is depends on the specific problem context.For this reason, in single-SNP models adjusting the p-values for FDR is important before selecting the final set of SNPs.While our proposed method is based on multi-SNP models, it may still need corrections to calibrate the potential of false discoveries.Finally, robustifying our proposed against data-level issues such as non-nuclear families, lack of individual-level data for some individuals warrant additional research.To this end, there is potential of adapting existing methods, such as Niu et al. [37], to the paradigm of e-values.Foundation (NSF) under grants 1737918, 1939916 and 1939956.

Figure 1 .
Figure 1.Density plots of bootstrap approximations for E * and E −j for all j in simulation setup, with s = 0.2, 0.3, 0.6, 1.

Table 1 .
(Top) Average True Positive (TP)/ True Negative (TN) rates for mBIC2, RFGLS+BH and the e-values method with E 1 and E 2 as evaluation maps and different values of t over 1000 replications, and (Bottom) Average Relaxed True Positive (RTP) and Relaxed True Negative (RTN) rates

Table 2 .
Table of analyzed genes and number of detected SNPs in them by the three methods

Table 3 .
Table of detected SNPs with known references