KW is a non-parametric statistical test used to assess whether there is a significant difference between the medians of two or more independent groups. Because the KW test is not originally designed to detect differences in variances, I have to convert the data into a measure of “dispersion” before using KW. Let e i j be the “covariate-adjusted residual” (obtained from Section 2.1) of individual j with genotype i, e i ∼ be the median of “covariate-adjusted residuals” within genotype group i. The deviation between individual j’s covariate-adjusted residual and his/her group median is D i j = e i j − e i ∼ (a measure of dispersion), where j = 1 , ⋯ , N and usually i = 1 , 2 , 3 (subscript j is sufficient to distinguish all N individuals; subscript i is supplementary to indicate the genotype group). Suppose r i j is the rank of D i j counting from all individuals across three genotype groups. The KW test statistic (Kruskal and Wallis, 1952) is listed in Equation 1 as follows, KW = N − 1 ∑ i = 1 M n i r i . ¯ − r ¯ 2 ∑ i = 1 M ∑ j = 1 n i r i j − r ¯ 2 (1) where N is the total sample size, n i is the number of individuals in genotype group i, M is the number of genotype groups (usually M = 3), r i . ¯ = ∑ j = 1 n i r i j / n i is the average of ranks in genotype group i, r ¯ = ∑ i = 1 M ∑ j = 1 n i r i j / N = 1 2 N + 1 is the average of ranks among all individuals. The KW test statistic (1) follows the chi-square distribution with the degrees of freedom M − 1 . All analyses in this work were conducted using R (version 4.3.1). The KW test was implemented with the R function “kruskal.test” in the R environment.

Although I have adjusted the covariates’ effects in Section 2.1, I still present the statistical model shown in QUAIL’s paper as follows (Miao et al., 2022). Through this model, people may get a complete understanding of Miao et al.’s methodology. Let the conditional quantile function (Q) of the trait Y at a SNP be Equation 2 as follows Q Y τ | g e n o t y p e = i , X = μ τ + i β τ + X α τ (2) where τ is the quantile ranging from 0 to 0.5, i = 1 , 2 , 3 , μ τ is the intercept, β τ is the coefficient of the quantile regression, X is the N × C matrix for C covariates, and α τ is the C-length vector for covariate effects. QUAIL measures β 1 − τ − β τ given τ ∈ 0 , 0.5 , i.e., the difference between the regression coefficients of the 1 − τ ^th and τ ^th quantiles. QUAIL tests the hypothesis H 0 : β Q I = 0 vs H 1 : β Q I ≠ 0 , where β Q I in Equation 3 is the quantile-integrated effect, i.e., β Q I = ∫ 0 0.5 β 1 − τ − β τ d τ (3) β Q I is the integral of β 1 − τ − β τ from τ = 0 to τ = 0.5 , which can be estimated by summing the areas of K rectangles with a length of β 1 − τ k ^ − β τ k ^ (where k = 1, … , K) and a width of 0.5 − 0 / K = 1 / 2 K . That is, β Q I can be estimated through β Q I ^ = ∑ k = 1 K β 1 − τ k ^ − β τ k ^ / 2 K . To have a more accurate estimation for β Q I , K is required to be a large number. In the QUAIL R code, Miao et al. used K = 100 as the default setting (Miao et al., 2022). The QUAIL R code was downloaded from GitHub (https://github.com/qlu-lab/QUAIL).

To search for SNPs that are associated with phenotypic variability, Marderstein et al. regressed D i j on the minor allele count at each SNP (i.e., 0, 1, or 2) (Marderstein et al., 2021). The regression model is D i j = α + β · i − 1 + ε j , where the predictor i − 1 is the minor allele count (i.e., 0, 1, or 2) and the random error term ε j follows a normal distribution with a mean of 0 and a variance of σ ε 2 . Because D i j is the deviation between individuals’ covariate-adjusted residuals and their group medians, this approach is called the deviation regression model (DRM). DRM tests whether the regression coefficient β is statistically significant, i.e., H 0 : β = 0 vs H 0 : β ≠ 0 . The DRM R code was downloaded from GitHub (https://github.com/drewmard/DRM).

The BF test statistic is listed in Equation 4 as follows, BF = ∑ i = 1 M n i D i . ¯ − D ¯ 2 / M − 1 ∑ i = 1 M ∑ j = 1 n i D i j − D i . ¯ 2 / N − M (4) where D i . ¯ = ∑ j = 1 n i D i j / n i is the average of deviations in genotype group i, D ¯ = ∑ i = 1 M ∑ j = 1 n i D i j / N is the average of deviations in all individuals. The BF test statistic (4) follows the F distribution with the numerator degrees of freedom M − 1 , and the denominator degrees of freedom N − M (Brown and Forsythe, 1974). The BF test is equivalent to the ANOVA F test for regressing D i j on two indicator variables that separate the three genotype groups of an SNP. Therefore, BF treats genotypes as a categorical scale (by using two indicator variables), while DRM regards genotypes as a continuous scale (by coding genotypes as 0, 1, or 2).

To implement the BF test, I used the R function “leveneTest” while specifying “center = median” in the “car” package (version 3.1–3) (Fox and Weisberg, 2018). The Levene’s and the BF tests are both used to test the homoscedasticity across different groups. The BF test is more robust to outliers than Levene’s test because it uses deviations from the group median, while Levene’s test uses deviations from the group mean.

DGLM fits two generalized linear models (GLM) that model genetic effects on the mean and dispersion of the covariate-adjusted residuals (Zhang and Bell, 2024). It regresses the covariate-adjusted residual e i j on the minor allele count at each SNP (i.e., 0, 1, or 2) while allowing the dispersion of e i j varies with the minor allele count. The regression model is e i j = α + β · i − 1 + ε j , where the random error term ε j follows a normal distribution with a mean of 0 and a variance of σ i 2 (allowing each genotype group has its variance). The DGLM tests whether the variance of e i j varies with the predictor (i.e., the minor allele count at each SNP). The DGLM method was implemented with the “dglm” R package (version: 1.8.6).

Figure 1 presents the QQ plots (the left column) and power under α = 5E-8 (columns 2–5) for five vQTL tests (N = 30,000; with SNP main effects). Supplementary Figure S1 demonstrates the results given N = 30,000, and SNP main effects do not exist. Supplementary Figure S2 (with SNP main effects) and S3 (without SNP main effects) show the results when N = 147,836. The rank-based inverse-normal transformation (INT) can convert data to a normal distribution. However, previous vQTL studies found that this INT transformation led to inflation in false positive rates (FPR) (Marderstein et al., 2021; Miao et al., 2022; Wang et al., 2019). The FPR was calculated by dividing the incorrectly classified negatives by the total negatives. To evaluate INT’s performance in our simulation setting, we equipped it with DGLM.

The statistical power was calculated by dividing the correctly classified positives by the total positives (each scenario was simulated 20,000 times). When the phenotype is normally distributed, DGLM and DGLM_INT were the most powerful methods (Figure 1; Supplementary Figure S1–S3C–E). Among the two non-parametric tests, QUAIL was superior to KW for normally distributed phenotypes (Figure 1; Supplementary Figure S1–S3C–E). When the phenotype is a kurtotic distribution, DGLM and BF had inflation in FPR (Figure 1; Supplementary Figure S1–S3F). KW became the most outstanding test among the four valid methods that preserve a suitable FPR (Figure 1; Supplementary Figure S1–S3F–J). When the phenotype is skewed, KW had a slight inflation in FPR, whereas DGLM and DGLM_INT had a severe inflation in FPR (Figure 1; Supplementary Figure S1–S3K). The other three valid methods (QUAIL, DRM, and BF) had similar power (Figure 1; Supplementary Figure S1–S3L–O).

By comparing the results of different simulation scenarios, I concluded three points: (1) When the MAF for SNP1 or SNP2 was larger, the percentage of the variance explained by the interaction effect increased, and the power of each method was enhanced. (2) When the sample size was enlarged from N = 30,000 to N = 147,836, the power performances of all methods were greatly improved. (3) Lastly, by comparing Figure 1 with Supplementary Figure S1, S2 with Supplementary Figure S3, the power of each method was higher in the presence of SNP main effects than in the absence of SNP main effects.

Table 1 shows the computation time (in seconds) for each simulation replicate. I measured the execution time in R (version 4.3.1) on a Linux system running at 3.6 GHz and 32 GB of RAM. Although QUAIL is a valid test with appropriate FPR under three trait distributions, it takes much more computation time than other methods. The other non-parametric test, KW, requires only 1/70–1/90 of QUAIL’s execution time. The computation time of the four parametric tests is also reasonable. The time of computation needed for these tests is ordered as DRM < BF < KW < DGLM ≈ DGLM_INT <<<<< QUAIL. Table 2 summarizes the performance of each method in the simulation on FPR, power, and computational time under different scenarios.

When conducting actual genome-wide data analysis, investigators had to exclude cryptic relatedness among individuals (Marees et al., 2018). The TWB investigated cryptic relatedness among participants with the KING (Kinship-based INference for GWAS) software (Manichaikul et al., 2010). I removed the person with higher missing genotype rates for each first- or second-degree relative pair. This step excluded 28,928 from the 147,836 TWB individuals, and 118,908 remained in the vQTL analysis.

I analyzed four lipid traits, including high-density lipoprotein cholesterol (HDL), low-density lipoprotein cholesterol (LDL), triglyceride (TG), and total cholesterol (TCHO). Supplementary Figure S16 shows the histograms of the four lipid traits. All four traits are right-skewed with positive skewness values, especially TG (skewness = 2.16). I first removed SNP’s main effect and covariates’ effects from each phenotype. Specifically, I extracted the residuals of regressing each trait on G 1 , j , G 2 , j (two indicator variables separating the three genotype groups of a SNP), and 17 covariates, including sex (male vs female), age (in years), BMI (in kg/m²), current smoking status (yes vs no), current drinking status (yes vs no), performing physical exercise (yes vs no), educational attainment (integer ranging from 1 to 7), and 10 ancestry PCs. Current smoking was defined as “having smoked cigarettes for at least 6 months when joining the TWB”. Drinking indicated “having a weekly intake of more than 150 mL of alcoholic beverages for at least 6 months when joining the TWB”. Regular exercise meant “performing exercise lasting for 30 min thrice a week”. Educational attainment was coded as an integer ranging from 1 to 7: 1 (illiterate), 2 (no formal education but literate), 3 (primary school graduate), 4 (junior high school graduate), 5 (senior high school graduate), 6 (college graduate), and 7 (Master’s or higher degree).

Sex, age, and BMI significantly influence people’s lipid profiles (Beyene et al., 2020). Smoking impairs lipid metabolism by decreasing LDL receptor expression (Ma et al., 2020). Alcohol consumption disturbs lipid metabolism by increasing adipose tissue lipolysis and causing ectopic fat deposition in the liver (Steiner and Lang, 2017). Regular physical activity has been found to improve lipoprotein-lipid profiles (Lin, 2021). Moreover, lower educational attainment is usually linked to worse lipid profiles (Lara and Amigo, 2018). Because these seven covariates are associated with lipid traits, I adjusted them in all vQTL methods. Furthermore, the first 10 ancestry PCs are adjusted in genetic analyses to avoid population stratification (Uffelmann et al., 2021).

Based on the simulation results, I chose the KW test as the primary vQTL method. This non-parametric approach is robust to outliers by transforming trait values into ranks (Figure 1; Supplementary Figure S1–S3F). Moreover, although KW has a slight inflation in FPR for skewed trait distributions (Figure 1; Supplementary Figure S1–S3K), investigators can perform a follow-up regression analysis to check whether GxE and GxG exist (i.e., the so-called “direct GxE or GxG analysis”). Furthermore, KW needs only 1/70–1/90 of the other non-parametric competitor’s (QUAIL) execution time. Because SNPs with MAFs <0.05 are difficult to replicate in GxG or GxE findings (Lin, 2024; Wang et al., 2019), I only analyzed 2,580,790 common SNPs with MAFs ≥ 0.05.

With the PLINK clumping procedure (Purcell et al., 2007), I detected 10 independent HDL-vQTLs (KW p-value < 5E-8) with linkage disequilibrium (LD) measure r ² < 0.01. Moreover, 6, 12, and 6 independent vQTLs were identified from LDL, TG, and TCHO, respectively. Four SNPs were found to be vQTLs of multiple lipid traits, including rs4704208 (near HMGCR), rs11748027 (in ANKDD1B), rs483082 (in APOC1), and rs662799 (in APOA5). Table 3 summarizes the 30 (=10 + 6+12+6–4) unique vQTLs of the four lipid traits.

From the regression (or variable selection) perspective, it is important to keep the hierarchical structure between main and interaction effects (Zhou et al., 2021). Following previous vQTL research (Wang et al., 2019), after identifying vQTLs, I performed a “direct GxE analysis” by regressing the trait on the minor allele count of each vQTL (0, 1, or 2), seven environmental factors (Es), and the product (interaction) term between the minor allele count and one of the seven Es. To remove population stratification, I also adjusted the top 10 ancestry PCs in this regression. The seven Es are listed as the horizontal axis of Figure 2, including SEX (female vs male), SPO (performing regular exercise, yes vs no), EDU (educational attainment, an integer from 1 to 7), AGE (chronological age, in years), BMI (in kg/m²), DRK (alcohol consumption, yes vs no), and SMK (cigarette smoking status, yes vs no). As Westerman et al. (Westerman et al., 2022) indicated, some vQTLs were “pleiotropic” concerning phenotypic variability. Table 2 also shows that four loci are vQTLs shared by multiple lipid traits. Instead of performing the direct GxE and GxG analysis for respective vQTLs of each lipid trait (i.e., 10 vQTLs for HDL, 6 vQTLs for LDL, 12 vQTLs for TG, and 6 vQTLs for TCHO), I analyzed all 30 vQTLs to have a more comprehensive picture of the four lipid traits.

Because TG had more vQTLs compared with the other three traits, I present the phylogenetic heat map of the direct GxE analysis for TG in Figure 2. The x-axis of Figure 2 lists the 7 Es. The magnitude of the value in Figure 2 represents–log10 (two-sided p-value of the SNP-E interaction), which is always positive. However, I deliberately added a positive/negative sign before the magnitude. A positive sign indicates that E exacerbates the vQTLs’ effects. In contrast, a negative sign suggests that the E attenuates the vQTLs’ effects. Females have more attenuated TG genetic effects than males (blue color in SEX column), whereas higher BMI, alcohol consumption, and cigarette smoking lead to more substantial TG genetic effects (red color in BMI, DRK, and SMK columns). Moreover, the phylogenetic heat maps for the remaining three traits are demonstrated in Supplementary Figures S17–S19.

I further performed a “direct GxG analysis” by regressing the trait on the minor allele counts (0, 1, or 2) of two SNPs, the seven Es mentioned above, and the product (interaction) term of the minor allele counts of the two SNPs. Similarly, I adjusted the top 10 ancestry PCs in this regression. Figure 3 shows the phylogenetic heat map of the direct GxG analysis for TG. The maps for the other three lipid traits can be found in Supplementary Figures S20–S22.

Regarding TG, the interaction between rs438811 (chr. 19, in APOC1) and rs662799 (chr. 11, in APOA5) is the most significant among the 30 2 = 435 GxG tests constructed by the 30 vQTLs (p = 1.3E-15). The magnitude of the value in Figure 3 represents–log10 (two-sided p-value of the SNP-SNP interaction), which is always positive. The most significant SNP-SNP interaction is between rs438811 (chr. 19, in APOC1) and rs662799 (chr. 11, in APOA5), where p = 1.3E-15. SNP rs483082 is highly correlated with rs438811 (r ² = 0.99), and rs3741297 is in obvious linkage disequilibrium with rs662799 (r ² = 0.19). Therefore, despite four black cells in Figure 3 indicating four significant SNP-SNP interactions, I only highlighted rs438811-rs662799 interaction. As analyzed by the direct GxG regression model, each T allele at rs438811 is associated with 10.8 mg/dL (p = 7.0E-117), and each G allele at rs662799 is associated with 26.4 mg/dL (p ∼ 0, extremely low p-value). Possible values of the interaction (product) term between rs438811 and rs662799 are 0, 1 (= 1 × 1 ), 2 (= 1 × 2 or 2 × 1 ), and 4 (= 2 × 2 ). Each point of the interaction term is associated with 5.9 mg/dL (p = 1.3E-15). This result indicates a synergistic interaction between rs438811 and rs662799: T allele at rs438811 and G allele at rs662799 working together exerts a more substantial impact on TG than they would have on their own.