Functional Regression Models for Epistasis Analysis of Multiple Quantitative Traits

To date, most genetic analyses of phenotypes have focused on analyzing single traits or analyzing each phenotype independently. However, joint epistasis analysis of multiple complementary traits will increase statistical power and improve our understanding of the complicated genetic structure of the complex diseases. Despite their importance in uncovering the genetic structure of complex traits, the statistical methods for identifying epistasis in multiple phenotypes remains fundamentally unexplored. To fill this gap, we formulate a test for interaction between two genes in multiple quantitative trait analysis as a multiple functional regression (MFRG) in which the genotype functions (genetic variant profiles) are defined as a function of the genomic position of the genetic variants. We use large-scale simulations to calculate Type I error rates for testing interaction between two genes with multiple phenotypes and to compare the power with multivariate pairwise interaction analysis and single trait interaction analysis by a single variate functional regression model. To further evaluate performance, the MFRG for epistasis analysis is applied to five phenotypes of exome sequence data from the NHLBI’s Exome Sequencing Project (ESP) to detect pleiotropic epistasis. A total of 267 pairs of genes that formed a genetic interaction network showed significant evidence of epistasis influencing five traits. The results demonstrate that the joint interaction analysis of multiple phenotypes has a much higher power to detect interaction than the interaction analysis of a single trait and may open a new direction to fully uncovering the genetic structure of multiple phenotypes.


Introduction
In the past several years, we have witnessed remarkable progresses in the development of methodologies for identification of epistasis that detect deviation from summation of genetic additive effects for a quantitative trait [1]. The classical approach to epistasis analysis is a single variant test. The epistasis is typically evaluated by testing interaction between a pair of variants one at a time. The classical methods for epistasis tests are originally designed to detect epistasis for common variants and are difficult applied to rare variants due to multiple testing problems and the low power to detect interaction. To overcome the critical barrier in interaction analysis for rare variants, instead of testing each pair of variants individually, group interaction tests that evaluate cumulative interaction effects of multiple genetic variants in a region or gene have recently been developed. Regression-based methods [2][3][4][5][6][7][8], haplotype-based methods [9][10][11][12][13][14][15], and machine learning-based methods [16][17][18][19][20] are proposed for epistasis analysis.
The classical statistical methods for interaction analysis have mainly tested association with single traits, one time analyzing one trait [21]. However, multiple phenotypes are highly correlated. More than 4.6% of the SNPs and 16.9% of the genes in previous genome-wide association studies (GWAS) are reported to be significantly associated with more than one trait [22]. These results demonstrate that genetic pleiotropic effects likely play a crucial role in the molecular basis of correlated phenotypes [23][24][25][26]. Joint epistasis analysis of multiple complementary traits will increase statistical power to unravel the interaction structure of multiple phenotypes [27,28]. Despite their importance in understanding genetic mechanism underlying the complex diseases, the statistical methods for identifying epistasis in multiple phenotypes have been less developed [1]. The interaction analyses for multiple phenotypes have been limited to common variants in carefully controlled experimental crosses [29,30]. Simultaneously analyzing interactions for multiple phenotypes in humans poses enormous challenges for methodologies and computations.
Purpose of this paper is to develop a general analytic framework and novel statistical methods for simultaneous epistasis analysis of multiple correlated phenotypes. To unify the approach to epistasis analysis for both common and rare variants, we take a genome region (or gene) as a basic unit of interaction analysis and use all the information that can be accessed to collectively test interaction between all possible pairs of SNPs within two genome regions (or genes). Functional data analysis is used to reduce the dimension of next-generation sequencing data. Specifically, genetic variant profiles that will recognize information contained in the physical location of the SNP are used as a major data form. The densely typed genetic variants in a genomic region for each individual are so close that these genetic variant profiles can be treated as observed data taken from curves [8,31]. Since standard multivariate statistical analyses often fail with functional data [32] we formulate a test for interaction between two genomic regions in multiple quantitative trait analysis as a multiple functional regression (MFRG) model [33] with scalar response. In the MFRG model the genotype functions (genetic variant profiles) are defined as a function of the genomic position of the genetic variants rather than a set of discrete genotype values and the quantitative trait is predicted by genotype functions with their interaction terms. By functional principal component analysis, the genotype functions are expanded as a few functional principal components (FPC) and the MFRG model is transformed to the classical multivariate regression model (MRG) in which FPC scores are taken as variates. Statistics are developed in this publication which can be applied to pairwise interaction tests and gene-based interaction tests for multiple phenotypes. By investigating SNP-SNP interactions or gene-gene interactions that are shared across multiple traits, pleiotropic epistasis can be studied.
To evaluate performance for multiple traits epistasis analysis, large scale simulations are used to calculate the Type I error rates of the MFRG for testing interaction between two genomic regions with multiple phenotypes and to compare power with multivariate pair-wise interaction analysis and single trait interaction analysis by functional regression (FRG) model. To further evaluate performance, the MFRG for epistasis analysis is applied to five traits: high density lipoprotein (HDL), low density lipoprotein (LDL), total cholesterol, systolic blood pressure (SBP), and diastolic blood pressure (DBP), from exome sequence data from the NHLBI's Exome Sequencing Project (ESP) to detect pleiotropic epistasis.

Methods
Assume that n individuals are sampled. Let y ik , k = 1,2,. . .,K, be the k-th trait values of the i-th individual. Consider two genomic regions [a 1 , b 1 ] and [a 2 , b 2 ]. Let x i (t) and x i (s) be genotypic functions of the i-th individual defined in the regions [a 1 , b 1 ] and [a 2 , b 2 ], respectively. Let y i = [y i1 ,. . .,y iK ] T be the vector of the trait values measured on the i-th individual. Let t and s be a genomic position in the first and second genomic regions, respectively. Define a genotype profile x i (t) of the i-th individual as where M and m are two alleles of the marker at the genomic position t. Recall that a regression model for interaction analysis with the k-th trait is defined as where μ k is an overall mean of the k-th trait, z kd is the coefficient associated with the covariate ν d , α kj is the main genetic additive effect of the j-th SNP in the first genomic region for the k-th trait, β kl is the main genetic additive effect of the l-th SNP in the second genomic region for the k-th trait, γ kjl is an additive × additive interaction effect between the j-th SNP in the first genomic region and the l-th SNP in the second genomic region for the k-th trait; x ij and z il are indicator variableS for the genotypes at the j-th SNP and the l-th SNP, respectively; ε ik , k = 1,..,K are independent and identically distributed normal variables with mean of zero and covariance matrix S. Similar to the multiple regression models for interaction analysis with multiple quantitative traits, the functional regression model for a quantitative trait can be defined as where α 0k is an overall mean, z kd is defined as before, α k (t) and β k (s) are genetic additive effects of two putative QTLs located at the genomic positions t and s, respectively; γ k (t,s) is the interaction effect between two putative QTLs located at the genomic positions t and s for the k-th trait, k = 1,. . .,K, x i (t) and x i (s) are genotype profiles, and ε ik are independent and identically distributed normal variables with mean of zero and covariance matrix S. Consider covariates in the model (2) allows incorporating PCA scores for population stratification, sex, age, BMI and other biomarkers into the model.

Estimation of Interaction Effects
We assume that both phenotypes and genotype profiles are centered. The genotype profiles x i (t) and x i (s) are expanded in terms of the orthonormal basis function as: x ij j ðtÞ and where ϕ j (t) and ψ l (s) are sequences of the orthonormal basis functions. The more number of variants in the genes the more accurate the eigenfunction expansion. If the number of variants is less than 3 the eigenfunction expansion of the genotypic profiles is impossible. MFRG can only be used for gene with more than 3 variants. In practice, numerical methods for the integral will be used to calculate the expansion coefficients. Substituting Eq (3) into Eq (2), we obtain (Appendix) The parameters α kj , β kl and γ kjl are referred to as genetic additive and additive × additive effect scores for the k-th trait. These scores can also be viewed as the expansion coefficients of the genetic effect functions with respect to orthonormal basis functions.
Then, Eq (4) can be approximated by (Appendix) where W ¼ ½ e ν x η Γ and B ¼ . Therefore, we transform the original functional regression interaction model into the classical multivariate regression interaction model by eigenfunction expansions. All methods for multivariate regression interaction analysis can directly be used for solving problem (5).
The standard least square estimators of B and the variance covariance matrix S are, respectively, given byB Denote the last JL row of the matrix (W T W) −1 W T by A. Then, the estimator of the parameter γ is given byγ The vector of the matrix γ can be written as By the assumption of the variance matrix of Y, we obtain the variance matrix of vec(Y): Thus, it follows from Eqs (9) and (10) that Test Statistics An essential problem in genetic interaction studies of the quantitative traits is to test the interaction between two genomic regions (or genes). Formally, we investigate the problem of testing the following hypothesis: which is equivalent to testing the hypothesis: Define the test statistic for testing the interaction between two genomic regions [a 1 , b 1 ] and [a 2 , b 2 ] with K quantitative traits as Then, under the null hypothesis H 0 : γ = 0, T I is asymptotically distributed as a central w 2 ðKJLÞ distribution if JL components are taken in the expansion Eq (3).
Group tests often make implicit homogeneity assumptions where all putatively functional variants within the same genomic region are assumed to have the same direction of effects. However, in practice, the variants with opposite directions of effects will be simultaneously presented in the same genomic region. MFRG can efficiently use information of both risk and protective variants and allow for sign and size heterogeneity of genetic variants. In general, the trait increasing and decreasing variants will be present in different locations in the genomic region. Information of trait increasing and decreasing variants usually will be reflected in different eigenfunctions and hence will be included in different functional principal component scores. The MFRG test statistic is essentially to summarize the square of the functional principal component scores. Therefore, the opposite effects of trait increasing and decreasing variants on the phenotype will not compromise each other in the MFRG test statistics. The MFRG statistics automatically take the opposite effects of the trait increasing and decreasing variants on the phenotype into account and do not require additional computations. MFRG will take the sign and size heterogeneity of the variants into account and be less sensitive to the presence of variants with opposite directions of effect.
We can also develop likelihood ratio-based statistics for testing interaction.
we can write the model as Under H 0 : γ = 0, we have the model: The estimators will bê The likelihood for the full model and reduced model are, respectively, given by nK=2 jŜj n=2 and nK=2 jŜ 1 j n=2 : The likelihood-ratio-based statistic for testing interaction between two genomic regions with multivariate traits is defined as Under the null hypothesis H 0 : γ = 0, T IΛ is asymptotically distributed as a central w 2 ðKJLÞ distribution if JL components are taken in the expansion Eq (3).

Simulation Model for Type 1 Error Rate Calculation
The genetic models for simulations to calculate Type 1 error rates of the tests are briefly given below. We first assume the model with no marginal effects for all traits: where Y i = [y i1 ,. . .,y ik ], μ = [μ 1, . . .,μ k ], and ε i is distributed as Then, we considered the model with marginal genetic effect (additive model) at one gene: where P j is a frequency of the allele A j , r k is a risk parameter of the k-th trait which was randomly selected from 1.1 to 1.6. The risk parameter affect the genetic effects and is used to control the contribution effort by genotype to the phenotype. The risk parameter influences the relative magnitude of the genetic effects. f 0 is a baseline penetrance and set to 1 and ε are defined as before.
Finally, we consider the model with marginal genetic effects (additive model) at both genes: P j and q l are frequencies of the alleles A j and B l , respectively, r pk and r qk are risk parameters of the k-th trait for the SNPs in the first and second genes, respectively, and randomly selected from 1.1 to 1.6, f 0 is a baseline penetrance and set to 1 and ε are defined as before.

Null Distribution of Test Statistics
To examine the null distribution of test statistics, we performed a series of simulation studies to compare their empirical levels with the nominal ones. We calculated the Type I error rates for rare alleles, and common alleles. To make simulations more close to real whole exome sequencing data, we generated 50,000 datasets consisting of 1,000,000 chromosomes randomly sampled from the NHLBI's Exome Sequencing Project (ESP) with 2,016 individuals and 18,587 genes. Each dataset included randomly selected a pair of genes from sequenced 18,587 genes. We randomly selected 20% of SNPs from each gene as causal variants. The number of sampled individuals from populations of 1,000,000 chromosomes ranged from 1,000 to 5,000. For each dataset, we repeated 5,000 simulations. We presented average type I error rates over 50,000 randomly selected pairs of genes from whole exome sequencing ESP dataset. Table 1 and S1 and S2 Tables summarized the average Type I error rates of the test statistics for testing the interaction between two genes with no marginal effect and consisting of only rare variants with 5 traits, 2 traits and 10 traits, respectively, over 50,000 pairs of genes at the nominal levels α = 0.05, α = 0.01 and α = 0.001. Table 2 and S3 and S4 Tables summarized the average Type I error rates of the test statistics for testing the interaction between two genes with marginal effect at one gene consisting of only rare variants with 5 traits, 2 traits and 10 traits, respectively, over 50,000 pairs of genes at the nominal levels α = 0.05, α = 0.01 and α = 0.001. Table 3 and S5 and S6 Tables summarized the average Type I error rates of the test statistics for testing the interaction between two genes with marginal effect at both genes consisting of only rare variants with 5 traits, 2 traits and 10 traits, respectively, over 50,000 pairs of genes at the nominal levels α = 0.05, α = 0.01 and α = 0.001. For common variants, we summarized the average Type I error rates of the test statistics for testing the interaction between two genes with marginal effect at both genes consisting of only common variants with 5 traits, 2 and 10 traits, respectively, over 10 pairs of genes at the nominal levels α = 0.05, α = 0.01 and α = 0.001, in Table 4 and S7 and S8 Tables, respectively. The statistics for testing interaction between two genomic regions with only common variants have the similar Type 1 error rates in the other two scenarios: with marginal genetic effects at one gene or without marginal genetic effects at two genes. These results clearly showed that the Type I error rates of the MFRG-based test statistics for testing interaction between two genes with multiple traits and common variants with or without marginal effects were not appreciably different from the nominal α levels. For the rare variants when the sample sizes increased to 5,000, the Type 1 error rates were still not appreciably different from the nominal levels.

Power Evaluation
To evaluate the performance of the MFRG models for interaction analysis of multiple traits, we used simulated data to estimate their power to detect interaction between two genes for two, four, five, six and ten quantitative traits. A true multiple quantitative genetic model is given as follows. Consider H pairs of quantitative trait loci (QTL) from two genes (genomic regions). Let Q h 1 and q h 1 be two alleles at the first QTL, and Q h 2 and q h 2 be two alleles at the second QTL, for the H pair of QTLs. Let u ijkl be the genotypes of the u-th individual with ij and g mu ijkl be its genotypic value for the m-th trait. The following multiple regression is used as a genetic model for the m-th quantitative trait: where g h mu ijkl is a genotypic value of the h-th pair of QTLs for the m-th quantitative trait and ε mu  Table). We assume that the genotypes at two loci affect a complex trait. Intuitively, Dominant OR Dominant model means that presence of risk allele at least one locus will cause the phenotype variation. Dominant AND Dominant model means that only when risk alleles at both loci are present the phenotype variation can be affected. Recessive OR recessive model indicates that when both risk alleles are at least present at one locus the phenotype variation can be observed. Threshold model implies that when two risk alleles at one locus and at least one risk allele at another locus are present, the phenotype variation will be observed. Recessive AND Recessive model is excluded due to low frequency of that condition with rare variants. The risk parameter r varies from 0 to 1.
We generated 2,000,000 chromosomes by resampling from 2,016 individuals of European origin with variants in random two genes selected from the NHLBI's Exome Sequencing Project (ESP). Two haplotypes were randomly sampled from the population and assigned to an individual. We randomly selected 20% of the variants as causal variants. A total of 2,000 individuals for the four interaction models were sampled from the populations. A total of 1,000 simulations were repeated for the power calculation.
The power of the proposed MFRG model is compared with the single trait functional regression (SFRG) model, the multi-trait pair-wise interaction test and the regression on principal components (PCs). For SNPs genotypes in each genomic region principal component analysis (PCA) were performed. The number of principal components for each individual which can explain 80% of the total genetic variation in the genomic region will be selected as the variables. Specifically, the principal component score of the i-th individual in the first and second genomic regions are denoted by x i1 ; . . . ; x ik 1 and z i1 ; . . . z ik 2 , respectively. The regression model for detection of interaction for the m-th trait is then given by x ij z il g mjl þ ε mi : Table 3. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effects at two genes consisting only rare variants with 5 traits over randomly selected 50,000 pairs of genes from the whole exome. The power of the MFRG is compared with the traditional point-wise interaction test which takes the following model: For a pair of genes, we assume that the first gene has k 1 SNPs, and the second gene has k 2 SNPs, then, the total number of all possible pairs is k = k 1 × k 2 . For each pair of SNPs, we calculated a statistic for testing pair-wise interaction T mjpair . Finally, the maximum of T mjpair : T max = max(T 1,1pair ,T 1,2pair ,. . .,T 1,kpair ,. . .,T M,1pair ,. . .,T M,kpair ) is computed. Only two genes include rare variants. These power curves are a function of the risk parameter at the significance level α = 0.05. Permutations in the point-wise interaction tests were used to adjust for multiple testing. In all cases, the two-trait FRG had the highest power to detect epistasis. We observed two remarkable features. First, two-trait test had higher power than the one-trait test. Second, the two-trait FRG had the highest power among all two-trait tests.      Dominant, Dominant AND Dominant, and Recessive OR Recessive, and r = 0.5 for Threshold models, respectively. Again, we observed that the power of the two-trait FRG was the highest.

Application to Real Data Examples
To further evaluate the performance, the MFRG for testing epistasis was applied to data from the NHLBI's ESP Project. Five phenotypes: HDL, LDL, total cholesterol, SBP and DBP were considered with a total of 2,016 individuals of European origin from 15 different cohorts in the ESP Project. No evidence of cohort-and/or phenotype-specific effects, or other systematic biases was found [34]. Exomes from related individuals were excluded from further analysis. We took the rank-based inverse normal transformation of the phenotypes [35] as trait values.
The total number of genes tested for interactions which included both common and rare variants was 18,587. The remaining annotated human genes which did not contain any SNPs in our dataset were excluded from the analysis. A P-value for declaring significant interaction after applying the Bonferroni correction for multiple tests was 2.89×10 −10 . Population stratification may inflate the test statistics. To reduce the inflation, the standard strategy is to adjust for population stratification via principal components. All the tests were adjusted for sex, age and population stratification via 5 principal components.
To examine the behavior of the MFRG, we plotted the QQ plot of the two-trait FRG test (Fig 7). The QQ plots showed that the false positive rate of the MFRG for detection of interaction in some degree is controlled.
A total of 91 pairs of genes which were derived from 85 genes showed significant evidence of epistasis with P-values < 2.7×10 −10 which were calculated using the MFRG model and simultaneously analyzing interaction of inverse normally transformed HDL and LDL (S10 Table). The top 30 pairs of significantly interacted genes with HDL and LDL were listed in Table 5. In Table 5 and S10 Table, P-values for testing interactions between genes by regression on PCA and the minimum of P-values for testing all possible pairs of SNPs between two genes using standard regression model simultaneously analyzed for the HDL and LDL and P-values for testing epistasis by the FRG separately against single trait HDL or LDL were also listed.
Several remarkable features from these results were observed. First, we observed that although pairs of genes showed no strong evidence of interactions influencing individual trait HDL or LDL, they indeed demonstrated significant interactions if interactions were simultaneously analyzed for correlated HDL and LDL. Second, the MFRG often had a much smaller P-value to detect interaction than regression on the PCA and the minimum of P-values of pairwise tests.
Third, pairs of SNPs between two genes jointly have significant interaction effects, but individually each pair of SNPs make mild contributions to the interaction effects as shown in Table 6. There were a total of 60 pairs of SNPs between genes CETP on chromosome 16 and GPR123 on chromosome 10 with P-values < 0.0488. None of the 60 pairs of SNPs showed strong evidence of interaction. However, a number of pairs of SNPs between genes CETP and GPR123 collectively demonstrated significant interaction influencing the traits HDL and LDL. Fourth, 91 pairs of interacting genes formed a network (Fig 8). The genes C5orf64 that had interactions with 19 genes, CSMD1 that had interactions with 20 genes, were hub genes in the network. 26 genes out of total 85 genes in the network were mainly located in 18 pathways. Each of 12 pathways included at least two interacting genes. However, the majority of interacting genes are located in different pathways. Among 18 pathways, calcium signaling pathway mediates the effect of LDL and plays a role in control of atherosclerosis susceptibility [36], LDL-cholesterol has multiple roles in regulating focal adhesion dynamics [37], LDL is involved in free radical induced apoptosis pathway [38], MAPK and JAK-STAT pathways are involved in dietary flavonoid protection against oxidized LDL [39], up-regulation of autophagy via AMPK/mTOR signaling pathway alleviates oxidized -LDL induced inflammation [40], PPARα holds a fundamental role in control of lipid homeostasis [41] and lectin-like ox-LDL receptor 1 mediates PKC-α/ERK/PPAR-γ/MMP pathway [42], HDL reduces the TGF-β1-induced collagen deposition [43], the Wnt pathway plays an important role in lipid storage and homeostasis [44], From the literatures, we found that both common and rare variants in CETP were associated with the HDL [45], CREBBP regulated LDL receptor transcription [46], PLTP was associated with HDL and LDL [47], TMEM57 was associated with serum lipid levels [48], SH2B3 was associated with LDL cholesterol [49]. It was also reported that CSMD1 was associated with multivariate phenotype defined as low levels of low density lipoprotein cholesterol (LDL-C < or = 100 mg/dl) and high levels of triglycerides (TG > or = 180 mg/dl) [50], associated with hypertension [51]. It was also reported that CSMD1 was associated with LDL and total cholesterol [52]. Next we analyzed five traits: HDL, LDL, SBP, DBP and TOTCHOL. Again, for each trait, inverse normal rank transformation was conducted to ensure that the normality assumption of the transformed trait variable was valid. To examine the behavior of the MFRG, we plotted QQ  plot of the test (S8 Fig). The QQ plots showed that the false positive rate of the MFRG for detection of interaction is controlled. A total of 267 pairs of genes which were derived from 160 genes showed significant evidence of epistasis influencing five traits with P-values < 1.96×10 −10 which were calculated using the MFRG model (S11 Table). Of them formed a largest connected subnetwork (Fig 9). The top 25 pairs of significantly interacted genes with five traits were listed in Table 7. We observed the same pattern as was observed for the two traits: HDL and LDL. 46 genes out of 160 genes in the networks were mainly located in 42 pathways including 15 signaling pathways. Among them, 14 pathways were in Fig 8. The interacting genes may be involved in the same biological pathway or in the different biological pathways. We observed 12 pathways, each of which contained at least two genes connected via interaction. However, the majority of interacting genes were not located in the same pathways.
Again, we observed that pairs of SNPs between two genes jointly have significant interaction effects, but individually each pair of SNPs might make mild contributions to the interaction effects as shown in S12 Table. There were a total of 6,766 pairs of SNPs between genes CSMD1 and FOXO1. S12 Table listed 101 pairs of SNPs with P-values < 0.049. The majority of the 101 pairs of SNPs showed no strong evidence of interaction. However, they collectively demonstrated significant interaction influencing five traits.
Among 42 pathways, in the previous sections we reported that 14 pathways were associated with HDL and LDL. From the literatures, we also know that unsaturated fatty acids stimulated the uptake of the LDL particles [53], PPAR signaling pathway was correlated with blood pressure [54], purine metabolism was associated with SBP [55], Wnt signaling pathway mediated cholesterol transportation [56], glycerolipid metabolism pathway was correlated with total cholesterol [57], focal adhesion pathway was involved in lipid modulation [58], Cell adhesion molecules was correlated with blood pressure [59].
We also observed from the literatures that a number of genes that appeared in the list of interacted genes with five traits had major genetic effects with single trait. Many reports showed that CETP, LIPC and LIPG were associated with HDL and LDL [60][61][62] and that MTHRR had known main effects for LDL [63] and blood pressure [64], NR1I3 for lipid metabolism [65], PLTP for LDL [66], [67], FOXO1 for LDL [68] and hypertension [69], SMAD9 for hypertension [70], and CSMD1 for SBP [51].

Discussion
Most genetic analyses of phenotypes have focused on analyzing single traits or, analyzing each phenotype independently. However, multiple phenotypes are highly correlated. Genetic variants can be associated with more than one trait. Genetic pleiotropic effects likely play a crucial role in the molecular basis of correlated phenotypes. To address these central themes and critical barriers in interaction analysis of multiple phenotypes, we shift the paradigm of interaction analysis from individual interaction analysis to pleiotropic interaction analysis and uncover the global organization of biological systems. MFRG was used to develop a novel statistical framework for joint interaction analysis of multiple correlated phenotypes. By large simulations and real data analysis the merits and limitations of the proposed new paradigm of joint interaction analysis of multiple phenotypes were demonstrated. The new approach fully uses all phenotype correlation information to jointly analyze interaction of multiple phenotypes. By large simulations and real data analysis, we showed that the proposed MFRG for joint interaction analysis of correlated multiple phenotypes substantially increased the power to detect interaction while keeping the Type 1 error rates of the test statistics under control. In real data analysis, we observed that although pairs of genes showed no strong evidence of interactions influencing individual trait, they indeed demonstrated significant interactions if interactions were simultaneously analyzed for correlated multiple traits.
Due to lack of power of the widely used statistics for testing interaction between loci and its computational intensity, exploration of genome-wide gene-gene interaction has been limited. Few significant interaction results have been observed. Many geneticists question the universe presence of significant gene-gene interaction. Our analysis showed that although the number of significantly interacted genes for single phenotype was small, the number of significantly interacted genes for multiple phenotypes substantially increased. Our results suggested that joint interaction analysis of multiple phenotypes should be advocated in future genetic studies of complex traits.
The interaction analysis for multiple phenotypes has been limited to common variants in carefully controlled experimental crosses and has mainly focused on the pair-wise interaction analysis. Although pair-wise interaction analysis is suitable for common variants, it is difficult to use to test interaction between rare and rare variants, and rare and common variants. There is an increasing need to develop statistics that can be used to test interactions among the entire allelic spectrum of variants for joint interaction analysis of multiple phenotypes. The MFRG utilizes the merits of taking genotype as functions and decomposes position varying genotype function into orthogonal eigenfunctions of genomic position. Only a few eigenfunctions that capture major information on genetic variation across the gene, are used to model the genetic variation. This substantially reduces the dimension in genetic variation of the data. The MFRG can efficiently test the interaction between rare and rare, rare and common, and common and common variants.
In both real data analysis of two phenotypes and five phenotypes, the interacted genes formed interaction networks. Hub genes in the interaction networks were also observed. These hub genes usually play an important biological role in causing phenotype variation.
An essential issue for interaction analysis of a large number of phenotypes is how to reduce dimension while fully exploiting complementary information in multiple phenotypes. The standard multivariate regression models for joint interaction analysis of multiple phenotypes do not explore the correlation structures of multiple phenotypes and reduce the dimensions of the phenotypes, and hence have limited power to detect pleotropic interaction effects due to large degrees of freedom. Data reduction techniques such as principal component analysis should be explored in the future interaction analysis of multiple phenotypes.
The results in this paper are preliminary. The current marginal approaches for interaction analysis cannot distinguish between direct and indirect interactions, which will decrease our power to unravel mechanisms underlying complex traits. To overcome these limitations, causal inference tools should be explored for the joint interaction analysis of multiple phenotypes. The purpose of this paper is to stimulate further discussions regarding great challenges we are facing in the interaction analysis of high dimensional phenotypic and genomic data produced by modern sensors and next-generation sequencing.  Table. Average type 1 error rates of the statistic for testing interaction between two genes with no marginal effect consisting of only rare variants with 2 traits over 10 pairs of genes. (DOCX) S2 Table. Average type 1 error rates of the statistic for testing interaction between two genes with no marginal effect consisting of only rare variants with 10 traits over 10 pairs of genes. (DOCX) S3 Table. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effect consisting of only rare variants with 2 traits over 10 pairs of genes. (DOCX) S4 Table. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effect at one gene consisting of only rare variants with 10 traits over 10 pairs of genes. (DOCX) S5 Table. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effects at two genes consisting of only rare variants with 2 traits over 10 pairs of genes. (DOCX) S6 Table. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effects at two genes consisting of only rare variants with 10 traits over 10 pairs of genes. (DOCX) S7 Table. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effects at two genes consisting of only common variants with 2 traits over 10 pairs of genes. (DOCX) S8 Table. Average type 1 error rates of the statistic for testing interaction between two genes with marginal effects at two genes consisting of only common variants with 10 traits over 10 pairs of genes. (DOCX) S9 Table. The interaction models: 0 and r stand for a quantitative trait mean given the genotypes.