Simulation Research on the Methods of Multi-Gene Region Association Analysis Based on a Functional Linear Model

Genome-wide association analysis is an important approach to identify genetic variants associated with complex traits. Complex traits are not only affected by single gene loci, but also by the interaction of multiple gene loci. Studies of association between gene regions and quantitative traits are of great significance in revealing the genetic mechanism of biological development. There have been a lot of studies on single-gene region association analysis, but the application of functional linear models in multi-gene region association analysis is still less. In this paper, a functional multi-gene region association analysis test method is proposed based on the functional linear model. From the three directions of common multi-gene region method, multi-gene region weighted method and multi-gene region loci weighted method, that test method is studied combined with computer simulation. The following conclusions are obtained through computer simulation: (a) The functional multi-gene region association analysis test method has higher power than the functional single gene region association analysis test method; (b) The functional multi-gene region weighted method performs better than the common functional multi-gene region method; (c) the functional multi-gene region loci weighted method is the best method for association analysis on three directions of the common multi-gene region method; (d) the performance of the Step method and Multi-gene region loci weighted Step for multi-gene regions is the best in general. Functional multi-gene region association analysis test method can theoretically provide a feasible method for the study of complex traits affected by multiple genes.


Introduction
Genetic analysis of rare variants is considered to be one of the most important components to compensate for the current deficiency of genetic variation, which has not yet been explained [1]. Although the lack of catalogs to speculate on the genotypes of rare variants and the high cost of sequencing technology have previously made it impossible to conduct very in-depth studies on rare variants [2,3], the development of high-throughput sequencing technology [4] has enabled scientists to obtain SNP data in a cheap and efficient way, as it contains a large amount of data on rare variants [5,6]. However, many previous tools and methods are designed for common variants, so there is still a lack of efficient and practical tools for rare variants association analysis. At present, single-marker association analysis is the most commonly used method of gene association analysis. However, if this method is directly applied to rare variants, it will be impossible to find loci with a moderate or low gene effect due to the limitations of single-marker association analysis [5,6]. The effect of a locus of a rare variant is small and not easily detected, and if the single-marker association analysis is used, many valuable association loci will be ignored. To better find the weak effect sites, some methods to make the weak effect more significant by concentrating the association information of the whole gene region were proposed: a) the method based on fold (Collapsing Methods) [7][8][9]. By directly compressing multiple loci into a new variable, the associated rare variants with weak effects distributed at multiple loci are aggregated to make it easier to find; b) the method based on kern [10][11][12]. When the variance of a set of random variables is 0, the set of variables is made up of the same value, and therefore the kernel by method only needs to check whether the variance component of the group of effect estimates corresponding to all genotype variables in the whole region is 0; c) the method based on functional data analysis [13][14][15][16]. Functional data analysis converts the discrete loci into a continuous variable through the basis function, and only the coefficients of the effect estimation function corresponding to the continuous variable need to be tested in association analysis. There is also a strategy to consider and examine multiple loci at the same time and determine the significance of each locus [17,18], which is more effective than single-marker gene association analysis because it considers the interrelationships between multiple loci. The simplest method is to use the multivariate linear model as the test model for the multi-gene locus test [17]. However, when only a small part of the multiple gene loci included in the test are related, the large degree of freedom of the uncorrelated loci will lead to the loss of power.
The above-mentioned methods, whether based on gene region information aggregation or multi-gene locus analysis, have their advantages and disadvantages. Moreover, through continuous improvement and innovation of experts and scholars, the shortcomings of these methods have been constantly overcome and their performance has become more and more excellent. Since the single-gene region approach can aggregate small effects, the multi-locus approach can improve the analysis power by considering the interrelationship between multiple variables. If we combine these two methods, we can expect to obtain a multi-region analysis method with the advantages of both methods. In addition, the actual situation of phenotypic tend to be controlled by a few gene regions. Some of these gene region effects are apparent, some are weak, and a strong effect is easy recognize. However, a weak effect can easily be concealed by a stronger effect, even considering that this part of the phenotypic is controlled by the effect of apparent genetic regions.
At present, some scholars have carried out research on the combination of the two ideas, i.e., the aggregation of genetic information in gene regions and the use of interrelationship among multi-gene regions: One of them is Turkmen and Lin [19], who further extended the statistical test PDT (Pedigree Disequilibrium Test [20]) and FBAT (Family-Based Association Test [21]), and proposed Block analysis methods. Firstly, the specific approach is to divide the gene sequence to be analyzed into a block-by-block in a certain way and assume that the variants within the region are interdependent, but that the relationship between the blocks is mutually independent; secondly, PDT or FBAT methods are used to analyze the loci of each block; thirdly, the results in the region are generalized by means of the squares sum of standardized variances. After the statistics of each small block are obtained, the squares sum of the statistics is calculated again. After the aggregation of the information twice, a statistic subject to chi-square distribution is obtained, which is used as the gene association analysis statistics of the large gene region composed of these blocks. This method assumes that the loci inside the block are interdependent and the information inside the block is aggregated by PDT and FBAT methods, assuming that the blocks are independent of each other. The other method is by Ayers and Cordell [22], who improved the two methods of Group Lasso and Group Sparse Lasso [23], enabling them to distinguish common and rare variants within a group. This method can test multiple groups at the same time, in which one group is treated as a variable and the relationship between different groups is considered.
The analysis method of a single gene region based on functional data is a method to express the high-density genetic markers as functional data through integration and analyze the region through a linear model. Many experts have shown that this is an effective way to improve the power of gene association [13,14,24]. If the functional linear model considers multiple gene regions at the same time, it can not only improve the power by considering the interaction between multi-gene regions, but also isolate the effects of each gene region through the characteristics of the linear model, making the gene regions with weak effects more obvious and easier to find. Therefore, we will explore the method of gene association analysis of multi-gene regions based on functional data analysis, hoping to find a method with higher power and better detection of association loci with weak effects, and provide some references for researchers interested in this field in the future.

Materials and Methods
Suppose that there are n individuals in a genetic population. The genome region [0, T] is constructed by SNP sequences t 1 ≤ t 2 · · · ≤ t M for genetic association analysis under the group structure and is not included. Let y i be the quantitative trait value of i-th individual and the population structure of the sample is not considered, so the traditional linear genetic model can be expressed as where µ 0 is the overall mean of the model, x ij is a genotype profile (if A and a represent a pair of alleles, then when the genotype of i-th individual is AA, x ij is taken as 2; when the genotype of i-th individual is Aa, it is taken as 1; when the genotype of i-th individual is aa, is 0). β j represents the effect coefficient of genetic marker, ε i ∼ N(0, σ 2 ), σ 2 is the environmental genetic variance, M is the number of genetic markers. With the increase of the number of genetic markers, the freedom degree of the model gradually increases, and the multicollinearity among variables becomes more and more serious, eventually leading to the reduction of estimation accuracy and power. This is especially true when the genetic markers are low-frequency variations. When the discrete variants are at ultrahigh density, the discrete variants in an interval are as continuous, and the functional linear model (FLM) can be used instead of the multiple linear genetic model: where ε i is an independent and normal distribution with zero mean and variance σ 2 , and [0, T] represents the genomic region under consideration. The discrete genetic markers x ij in Equation (1) are converted into continuous genetic markers function X i (t) in Equation (2). At this time X i (t) is a random function, and the effects of genetic markers β j are also converted into a continuous genetic effect function β(t). Step The functional linear genetic model of single-gene region is generally in the following form When the single-gene region is extended to the multi-gene region, the original linear genetic model becomes the following form for p-th (p = 1, 2, · · · , P) gene region. Every gene region is [0, T]. For every region, the lower bound of the interval is converted to zero, and the upper bound of the interval is converted to T. According to the method of functional data analysis, a set of basis functions ϕ p1 (t), ϕ p2 (t), · · · ϕ pK G (t) and the coefficient d pi1 , d pi2 , · · · , d piK G can be used to expand X pi (t) as Similarly, a set of basis functions φ p1 (t), φ p2 (t), · · · , φ pK β (t) and coefficient b p1 , b p2 , · · · , b pK β can be used to expand β p (t) as the expansion of X pi (t) and β p (t), the model (4) becomes Let , as well as (7) becomes where ϕ p1 (t), ϕ p2 (t), · · · ϕ pK G (t) and φ p1 (t), φ p2 (t), · · · , φ pK β (t) are a set of orthonormal basis. Usually, we choose the same basis function for ϕ p1 (t), ϕ p2 (t), · · · ϕ pK G (t) and φ p1 (t), φ p2 (t), · · · , φ pK β (t). Therefore, the above genetic model can be further simplified as Model (8) becomes model (9). At this point, the genetic model is transformed into the ordinary multiple linear regression model of Equation (9), for which variables can be screened by stepwise regression [25][26][27]. Because only the interrelationship between whole gene regions and traits is discussed, d pi represents the genetic information of p-th gene region, so d pi , which represents the whole information of a gene region, is considered to be added to the model as a "variable". p gene regions should be screened as p variables.
There are three ways of screening variables for regression: forward selection, backward selection, forward selection, and backward selection. Here, the backward selection method is performed, where all gene regions to be analyzed are put into the model at the beginning, and then some gene regions are removed step-by-step until a reduced model is obtained. AIC (Akaike Information Criterion) information criteria will be used as the basis for each step to determine which gene regions need to be removed from the model where Rss represents the squares sum of the residuals for the current model, and K represents the number of unknown variables, that is, the sum of the number of elements of all d pi . After deciding which d pi to remove, we hypothetically delete each d pi existing in the current model and calculate the AIC, which is made up of the rest of the gene region. We find the model corresponding to the minimum AIC and then proceed to the next step. By repeating the above steps, until deleting any of the gene regions in the model does not make the AIC of the model smaller, we now have a very reduced model. Finally, the partial F test commonly used in the multiple linear regression model was used to test each d pi in the model, and the corresponding p-value was calculated as the evaluation basis for the association between gene regions and quantitative traits.

Multi-SLoS
Lin et al. [28] proposed a locally sparse functional linear model (SLoS method, Smooth and Locally Sparse method). By adding fSCAD (Functional Smoothly Clipped Absolute Deviation) penalty function on the basis of smoothing penalty term, the functional linear model has the ability to identify the sparse part of the estimated effect valueβ(t) and compress the estimated value to null. In the paper, a single region is taken as an example, but a model for multi-regions was also proposed: According to the description of the paper, to estimate the corresponding β(t) = (β 1 (t), β 2 (t), · · · β P (t)) T , we only need to solve the corresponding loss function: Specific algorithms can be found in Lin et al. [28]. In order to distinguish the SLoS method for a single region, we refer to it as the Multi-SLoS method. The Multi-SLoS method has the ability of local sparse, that is, it can identify null and no-null in β(t). We hope to use the local sparsity ability of this model for gene association test, which involves the problem of the model test. That is, the significance test problem of individual gene regions in the model (2) is presented based on the Multi-SLoS method. The following is a detailed description of how to conduct the test based on the work of Lin et al. [28] for multi-gene region association analysis. The functional linear model of a single gene region can be transformed into a multivariate linear model of several variables (the number of basis functions plus intercept), and then directly test whether the estimated coefficients of the basis functions are all null [13,29,30]. Multi-regions association analysis can be done similarly. However, when multi-gene regions are tested, not only must more variables be tested, but also the degrees of freedom of the model should be adjusted. In addition, we found that the results obtained after adjusting the degrees of freedom of the Multi-SLoS method as follows would be more consistent with the features of the method and the actual results in the follow-up study of polygenic regions. The reasons for the adjustment of degrees of freedom are given below.
Let the null gene regions denote gene regions where estimated effect valuesβ p (t)(p = 1, 2, · · · , P) are all null, and non-null gene regions denote gene regions where estimated effect valueβ p (t) are not all null. The Multi-SLoS method will directly compress the estimated effect valuesβ p (t) of null gene regions to null and identify the non-null and null gene regions from p gene regions, regardless of whether the Multi-SLoS method is correct in distinguishing the non-null and null gene regions (the results can be seen in the later simulation). The regions whereβ p (t) is compressed to null have no effect on the estimated results, which means that these regions have been directly identified asβ p (t) = 0 at a certain estimation stage of the model and have no effect on the estimated model. Then, the degrees of freedom of these gene regions should be removed from the calculation of the model. Similarly, there are some sub-regions where the effect is also compressed to null in non-null gene regions, which means that these sub-regions also have no effect on the estimated model, and the degrees of freedom of these sub-regions should be deducted.
Combined with the adjustment of degrees of freedom, the partial F test of the Multi-SLoS method is given below.
Let B represent the set of subscripts of gene regions with non-null effect, and b represent one element of the set B. For the above linear genetic model, the following method is used to test the association between gene region b and quantitative traits.
(a) Calculate the sum of square residuals of the full model including all non-null gene regions:ŷ (b) Calculate the sum of square residuals of the reduced model excluding gene region a: (c) Adjustable degrees of freedom: The adjusted freedom degree of SSE (full) (freedom adj (full)) should be: the number of individuals (n)-the sum of the number of non-null basis function coefficients in all non-null gene regions-1.
The adjusted freedom degree of SSE (reduced)-SSE (full) (freedom adj (reduced)) should be: the number of non-null basis function coefficients in gene region a. ∼ F( f reedom adj (reduced), f reedom adj ( f ull)

Multi-Gene Region Weighted Method
In common gene association analysis, common variants can be easily identified if the associated loci contain both rare and common variants, but rare variants are difficult to detect because of their micro effects. For the association analysis of multi-gene regions, a similar situation is likely to occur-only the associated regions with rare variants are difficult to identify if the associated loci exist in the regions only with common variants and the regions only with rare variants at the same time. The common solution to this problem in gene-association analysis is to assign different weights to different types of variants. The same approach is used to assign different weights to different gene regions, by assigning weights to different types of gene regions to eliminate differences due to different allele frequencies rather than different degrees of association with phenotypic values.

Weighted SLoS (W-SLoS)
The SLoS approach has been applied to polygenic regions in Section 2.1.1, which asks the question: would it further improve the power if different weights were given to different types of gene regions? The loss function of the SLoS method is: It can be seen from the loss function that different gene regions can be assigned different weights by adjusting parameters γ p and p λ p . Therefore, based on the research of Lin et al. [28], we have made appropriate adjustments to the code in the slos package of R language provided by Lin et al., so that the method is not only theoretically feasible but also runs smoothly in the actual program. Finally, the Weighted SLoS method only increases the weight compared to the Multi-SLoS method, and the same statistical test can be used to test the significance of each gene region.

Multi-Gene Region Loci Weighted Method
Although it is possible to distinguish rare variants from common variants and then divide them into rare variants regions and common variants regions for analysis in the multi-gene region analysis, it is more common in the actual situation that both common variants and rare variants exist in a gene region to be analyzed. Therefore, the multi-gene region loci weighted method is proposed, which is a more general method of combining functional data analysis by assigning different weights to each locus within each region rather than to the gene region, as in Section 2.1.2.

Multi-Gene Region Loci Weighted
Step (LW-Step) Similar to Section 2.1.1, there are n individuals in a genetic population. The genome region [0, T] is constructed by SNPs sequences t 1 < t 2 < · · · < t M for genetic association analysis under no group structure. Accordingly, the genetic markers are x i1 , x i2 , · · · , x iM (i = 1, 2, · · · , n). Let y i be the quantitative trait value of the i-th individual and the population structure of the sample is not considered, and the traditional linear genetic model can be expressed as With the increase of the number of genetic markers, the functional linear model (FLM) can be used instead of the multiple linear genetic model A set of basis functions Genes 2022, 13, 455 8 of 22 According to functional analysis method [31], let X i = [x i1 , x i2 , · · · , x iM ] represent the gene data vector for the i-th individual, then, . . .
According to the functional data analysis methods, there are The coefficient d i is solved in a smooth way Here, R 2 is a penalty matrix, The In addition, where functional linear model the following can be obtained where Next, as in Belonogova et al. [32], a M × M diagonal matrix Θ was designed, where each element on the diagonal of the matrix corresponds to the weight of genotype data where a 1 , a 2 are the preset parameters, and MAF ij represents the j-th genotype frequency of the i-th individual. The diagonal matrix Θ is embedded to the simplified functional linear equation [32] y This gives us a weighted functional linear function. This is the method of a single gene region, corresponding to the functional linear equation of multiple gene regions, which can be All the assumptions about gene regions are the same before. According to the specific situation of each region, add a weight matrix Θ p to assign different weights to the loci in the region, then The above statement is used to better explain the method in theory, in fact, the actual processing is not so complicated as in theory. Returning to the fitting of genotype data, the problem of loci weighting can be viewed from another perspective. First, let X * i = X i Θ, and then expand X * i with the same functional smoothing parameters, so that These are the same smooth parameters, basis functions, number of basis functions, and nodes.
is the same as the expansion of X i (t) and X * i (t). The difference between X i (t) and X * i (t) is the weight matrix Θ. This result can be used to deduce the single gene loci weighted functional linear model as follows n. It can be seen from the derived results that the single gene loci weighted linear functional linear model can be understood as the weighted transformation of the original genotype data into new functional data X * i (t), and then the functional linear model can be established by using X * i (t). We extend this result into the multi-gene region weighted functional linear model, and all the assumptions are the same as the above section. The model becomes the ordinary multi-gene region functional linear model Then, And For the same reasons as in Step method, when ϕ p1 (t), ϕ p2 (t), · · · ϕ pK G (t) and φ p1 (t), φ p2 (t), · · · , φ pK β (t) are a set of orthonormal basis and the same basis function, the difference between different individual data is mainly reflected in the coefficients. Therefore, it is reasonable to simply take the coefficients of the basis function as a new variable as the basis of subsequent method operations, and simplify the model as follows The model becomes a multivariate linear model with the coefficients of the basis functions in the region as the new variables. The difference between the two methods compared with the Step method is that, at this time, the coefficients are obtained by functional data analysis of genotype data after weighting. The simplified model treats each gene region as a "variable" for stepwise regression and conducts a partial F test on the final simplified model to obtain the significance level of each gene region. The method is called the Multi-gene region loci weighted Step (LW-Step).

Multi-Gene Region Loci Weighted SLoS Method (LW-SLoS)
For the multi-gene region, the model becomes the ordinary multi-gene region functional linear model The SLoS method can be used to solve the model, and the statistical test method proposed in Multi-SLoS method can be used to test the model. Finally, in order to distinguish between other types of SLoS methods, we will call this method LW-SLoS, which means Loci Weight SLoS. In power simulation, the associated loci should be assumed as the target of the analysis and the quantitative traits should be simulated as the analysis objects. Therefore, in each simulation, five of the 25 gene regions splicing "multi-gene regions" are selected as the associated gene regions, and then three loci are randomly selected from each associated gene region as the associated loci. The generation of simulated traits adopted the additive effect model. At the same time, three different scenarios are made for the effect of the associated loci: Scenario I, all effect directions are positive; Scenario II, the effect direction of all associated loci are negative in two of the five associated regions; Scenario III, the effect direction of one locus in each associated region is negative. The absolute value of effect value is determined by the following effect model where MAF t j represents the minimum allele frequency of the t j -th genotype as the associated locus in the associated gene regions. In the false-positive proportional simulation, random numbers were generated through normal distribution N(0, 0.1) as the phenotypic values of the simulation, since it was assumed that no associated gene loci existed in the gene region.
The multi-gene regions that are composed of 25 gene regions are analyzed for every simulation. Each gene region is simulated 100 times under different association effect hypotheses. There are 2500 (100 × 25) gene regions analyzed. In each scenario, 5 gene regions are assumed to be associated regions, and 20 gene regions are assumed to be unrelated regions. That is, there are 500 associated gene regions and 2000 unrelated gene regions for every case. Under given significance level α based on every model, every method and every scenario, the number of significant gene regions is n 1 . The number of significant gene regions, but no significant gene regions, in fact, is n 2 . The power is n 1 500 and the false positive rates (Type I error rates) are n 2 2000 . In order to compare the Step and Multi-SLoS method of multi-gene region analysis with the single-gene region functional method, the SLoS method and FLM method are also performed in this simulation. The single-gene region method analyzes the sub-regions of the simulated multi-gene region one by one and then summarizes the results to test the multigene regions. The FLM method is proposed by Svishcheva et al. [14] as a functional gene region analysis method for family and population gene data. The author provides a package, "FREGAT", written in the R language, that contains the computer program for the method. Moreover, the method can also be used for genetic association analysis in populations with no family relationship. For SLoS method and Multi-SLoS method, the compression parameter and smoothing parameter are set as 0.1 and 0.1 in the common variants multigene regions simulation; as 0.01 and 0.01 in the rare variants multi-gene regions simulation; and as 0.05 and 0.05 for the hybrid variants multi-gene regions simulation and the hybrid variants multi-gene regions simulation. The fitting of SLoS and Multi-SLoS models requires the "SLoS" package of R language [28], but the statistical part needs to be supplemented by the R language code written by us. In the process of simulation, genotype data needs to be converted into functional data, all gene regions are smoothed by 25 Fourier basis functions, the number of nodes is equal to the number of variants in the gene regions, and the distance of nodes is equidistant. The number of basis functions of the effect function uses the default settings for the appropriate software. The basis functions and node settings are the same in subsequent computational simulations, except for additional instructions.

Simulation of Multi-Gene Region Weighted Method
We designed three different multi-gene regions: (a) rare variants multi-gene regions, i.e., all variants within the regions are rare variants (MAF < 0.01); (b) common variants multi-gene regions, i.e., all variants within the regions are common variants; (c) the hybrid variants multi-gene regions: 15 are common variants gene regions and 10 are rare variants gene regions. Here, a simulation of the multi-gene region weighted method is only aimed at hybrid variants multi-gene regions analyzed in the simulation of common multi-gene region method. Besides, there are two common variant gene regions and three rare variant gene regions for five associated regions. In addition to counting the power and false positive rates of common and rare variant gene sub-regions in multi-genic regions, the rest of the settings are similar with those simulations of the common multi-gene region method, in order to count the power and false positive rates of the common variants gene regions and the rare variants gene regions in the multi-gene regions respectively. We illustrate the rare variants gene regions as an example: the power of rare variants gene sub-regions is the total number of significant and true associated gene regions in the rare variants gene sub-regions divided by the total number of rare variants gene sub-regions in the multi-gene regions. False positive rates of rare variants gene sub-regions: the total number of significant but unrelated gene regions in rare variants gene sub-regions is divided by the total number of rare variants gene sub-regions in multi-gene regions.
Three methods are used in this simulation: Step, FLM, and Weighted SLoS (W-SLoS). During the simulation of the W-SLoS method, the smoothing parameters of common variants gene sub-regions are 0.02 and the compression parameters are 0.05. The smoothing parameter and compression parameter of rare variants gene sub-regions are 0.01 and 0.0025, respectively. All settings for the nodes and the effect of genetic variations are similar as the computer simulations for the common multi-gene region method (Section 2.2.1).

Simulation of Multi-Gene Region Loci Weighted Method
The multi-gene regions generated by the simulation in the multi-gene regions loci weighted analysis is composed of 10 mixed variants gene regions. The proportion of rare variants in the 10 regions of mixed variants are as follows: [0.7,0.8,0.6,0.95,0.9,0.95,0.7,0.9,0.8,0.6], respectively. Common variants and rare variants of every sub-region are random distributions for simulated multi-gene regions. The gene loci are weighted in the sub-regions. The rare variants and common variants in each gene region are generated in the same way as the common multi-gene regions simulations. This multi-gene region will be analyzed 100 times during the simulations. During the power simulation, the associated gene regions are preset as 2-th, 4-th, 5-th, 7-th, and 10-th gene regions, the proportion of rare variants in the corresponding gene regions is [0.8,0.95,0.9,0.7,0.6], and the 4-th and 5-th gene regions are adjacent. Although the associated gene regions are determined, the associated loci in each gene region are randomly extracted from loci with a minimum allele frequency of less than 0.02 within the corresponding gene regions in each simulation. In addition, the model and method used to simulate phenotypic values and the effect value scenario settings of the associated loci are the same as those in the common multi-gene regions simulations. Suppose that there are three scenarios for the effect of the associated loci: Scenario I, all effect directions of all associated loci are positive for gene loci; Scenario II, the effect directions of all associated loci are negative for the 4-th and 7-th gene region; Scenario III, choose a locus at random form associated loci and the effect direction of that locus is negative in each associated region.
For power and false positive rates in the simulation, it is assumed that in the set of regions that are significant under a certain condition, the number of regions that are truly correlated and significant is n 1 , while the number of uncorrelated regions and identified as significant is n 2 , and then the total power in the simulation is: n 1 500 , and the total false positive rates (Type I error rates) are: n 2 500 . Suppose that the number of regions identified as significant in the 2-th, 4-th, 5-th, 7-th, 10-th gene regions is n i1 , (i = 2, 4, 5, 7, 10), respectively. Then, the powers of the corresponding gene sub-regions are n i1 100 , (i = 2, 4, 5, 7, 10). Suppose that the number of regions in the 1-th, 3-th, 6-th, 8-th, 9-th gene regions identified as significant is: n i2 , (i = 1, 3, 6, 8, 9) (in fact, those gene regions are no significant), then the false positive rates of corresponding gene sub-regions are n i2 100 ,(i = 1, 3, 6, 8,9). Four methods are used in the simulation: Step, Loci Weighted Step (LW-Step), Multi-SLoS, and Loci Weighted SLoS (LW-SLoS). The smoothing and compression parameters of the Multi-SLoS and LW-Step methods are set to 0.001. The multi-gene region loci weighted method needs to set the weights of different loci. In this paper, it is realized by setting the weight matrix. Here, the weight matrix is set as the diagonal matrix, and the weight of the corresponding gene loci of the diagonal elements on the matrix is generated through the beta distribution: Beta(MAF i , 1, 10).

Results of the Common Multi-Gene Region Method
It can be seen from Table 1 that the power of the Step method is the highest among the four methods at each significance level in the rare variants multi-gene regions, and the power of the Multi-SLoS method is not higher than that of FLM and Step, but it is also slightly higher than that of SLoS; in hybrid variants multi-gene regions I and hybrid variants multi-gene regions II the power of the Step method is higher than the FLM method, but the power of the Multi-SLoS method is not higher than the SLoS method. Therefore, in terms of power performance, multi-gene region analysis has a comparative advantage in rare variants multi-gene regions, and the power of the Step method is the best in four simulated gene regions. By comparing the power on different types of gene regions, the SLoS and Multi-SLoS method have higher power for multi-gene regions with common variants; the Step method makes it less powerful in multi-gene regions that contain common variants; the FLM method is more effective in multi-gene regions consisting of only rare or common variants, but the performance of hybrid variants multi-gene regions I is very strange. Table 1. Simulation results of the power for four types of multi-gene regions regarding the common multi-gene region method.

Multi-Gene Region Multi-Gene Region Multi-Gene Region I Multi-Gene Region II
Step

Multi-SLoS FLM SLoS
Step Combined with the simulation results of false positive rates in Table 2, the results are further analyzed. In rare variants multi-gene regions, the Multi-SLoS method has a higher false positive rate than the SLoS method, and the Step method has a higher false positive rate than the FLM method under different effect directions and higher significance level. The false positive rates in both the Step method and Multi-SLoS method are very low in the common variants multi-gene regions. In the hybrid variants multi-gene regions II, the Step method and Multi-SLoS method compared to the FLM and SLoS methods have lower false positive rates. Therefore, the addition of common variants to the sub-regions of the multi-gene regions has the effect of reducing the false positive rates.
On the one hand, the multi-gene region Step method based on functional data analysis has the best power and performance in false positive rates. Indeed, it can better find some associated regions that cannot be found in single gene regions, especially some associated gene regions with micro effects. On the other hand, the Multi-SLoS method has no significant advantages over the SLoS method and needs further improvement and adjustment. Note: Scenario I-all effect directions are positive: Scenario II-the effect direction of all associated loci are negative in two of the five associated regions; Scenario III-the effect direction of one locus in each associated region is negative.

Results of Multi-Gene Region Weighted Method
As can be seen from the power simulation results in Tables 1 and 3, the highest power of the Multi-SLoS method is 76.2% and of the SLoS method is 78.4%. The lowest power of W-SLoS is 84% for the hybrid variant multi-gene region I (see Table 1). Therefore, the SLoS method with weighted multi-gene regions has a significant improvement in power. The FLM method still does not perform well in this case. Combined with the previous simulation of the FLM method, it has a good performance when the gene sub-region in the multi-gene regions is of the same type, but the performance is not good when the multi-gene regions are mixed with multiple types of gene regions. It can be known that this method may be more suitable for the detection of multi-gene regions with the same type of gene sub-regions. The Step method is still the most powerful of the four methods, both in terms of overall power and in gene regions with common or rare variants. Almost all of the common variant gene regions are identified, and rare variants gene regions have lower power than that of the common variants gene regions. In general, the three methods can more easily identify rare variants gene regions, but common variants gene regions still require higher power. Table 4 shows the simulation results of false positive rates. It is obvious that the false positive rates of FLM are larger than that of other methods. Among three methods, the false positive rates of the W-SLoS method are the lowest, followed by the Step method. In general, the performance of W-SLoS is better than that of Step on false positive rates from Table 4. It means the W-SLoS method is not as effective as the Step method for detecting gene regions, but it is more reliable for gene regions selected by the W-SLoS method.
By weighting the sub-regions, the performance of the multi-gene region SLoS method exceeds that of the single-gene region SLoS method. This indicates that weighting can further improve the power of the polygenic region analysis model. Table 3. Simulation results of the power for three types of multi-gene regions regarding the multi-gene region weighted method.

Common Variant Multi-Gene Region Rare Variant Multi-Gene Region Hybrid Variant Multi-Gene Region I W-SLoS
Step FLM W-SLoS Step FLM W-SLoS Step FLM Note: Scenario I-all effect directions are positive: Scenario II-the effect direction of all associated loci are negative in two of the five associated regions; Scenario III-the effect direction of one locus in each associated region is negative. Table 4. Simulation results of the false positive rate for three types of multi-gene regions regarding the multi-gene region weighted method.

Common Variant Multi-Gene Region Rare Variant Multi-Gene Region Hybrid Variant Multi-Gene Region I W-SLoS
Step FLM W-SLoS Step FLM W-SLoS Step FLM Note: Scenario I-all effect directions are positive: Scenario II-the effect direction of all associated loci are negative in two of the five associated regions; Scenario III-the effect direction of one locus in each associated region is negative. Tables 5 and 6 show the power and the false positive rates simulation results of the multi-gene region loci weighted method. In the power simulation results, there are three unweighted methods: FLM, Multi-SLoS and Step, and the highest power of these methods is 94.2%. The highest power of the two loci-weighted methods is 98.2%. This suggests that multi-gene regions loci weighted methods can indeed improve the power in this simulation. For false positive rates simulation, the false positive rates of LW-Step are higher than that of Step, but lower than that of the other three methods. The LW-SLoS method has much higher false positive rates than the LW-Step method, but the false positive rates of LW-SLoS approach that of the LW-Step method as the significance level increased gradually.

Results of Multi-Gene Region Loci Weighted Method
Since the multi-gene regions in each simulation are fixed, so are the associated gene sub-regions as well, and the power and false positive rates of each gene sub-region are calculated. Figures 1 and 2, respectively, show the power of the 2-th, 4-th, 5-th, 7-th, and 10-th associated gene regions, and the false positive rates of the 1-th, 3-th, 6-th, 8-th, and 9-th unassociated gene regions. Figure 1 shows that the power of different gene regions is different under different effect hypotheses, the proportion of common and rare variants in a gene region affects power of the region. Figure 2 shows that in general, the LW-SLoS method has higher false positive rates in the 1-th and 3-th gene regions. The proportions of rare variants in these two gene regions are 0.7 and 0.6, which are the second and the smallest in the five gene regions.
In general, the multi-gene region loci weighted method has an advantage in the multigene regions, where the proportions of various variants in each gene region are different. Moreover, the simulation results of power and false positive rates show that such results are not simply lowering the threshold, but it improves the power of the analysis at reasonable false positive rates. The simulation further analyses the power and false positive rates of different gene sub-regions in the same multi-gene regions. The results show that the presence of some common variants in the sub-gene regions could improve the power of the method and reduce false positive rates. in a gene region affects power of the region. Figure 2 shows that in general, the LW-SLoS method has higher false positive rates in the 1-th and 3-th gene regions. The proportions of rare variants in these two gene regions are 0.7 and 0.6, which are the second and the smallest in the five gene regions. In general, the multi-gene region loci weighted method has an advantage in the multi-gene regions, where the proportions of various variants in each gene region are different. Moreover, the simulation results of power and false positive rates show that such results are not simply lowering the threshold, but it improves the power of the analysis at reasonable false positive rates. The simulation further analyses the power and false positive rates of different gene sub-regions in the same multi-gene regions. The results show that the presence of some common variants in the sub-gene regions could improve the power of the method and reduce false positive rates.   Note: Scenario I-all effect directions are positive: Scenario II-the effect direction of all associated loci are negative in two of the five associated regions; Scenario III-the effect direction of one locus in each associated region is negative.

Discussion
In this paper, a total of five analysis methods are proposed for the analysis of multigene regions, which can be divided into three categories: the common multi-gene region method, the multi-gene region weighted method, and the multi-gene region loci weighted method. the Step method merged the two ideas together with the gene information of the region and the relationship among the gene regions. The simulation results showed that the power of the Step method is higher than that of the FLM and SLoS method for a single gene region, even for that of the region-weighted SLoS method. It means the power is improved for considering the relationship among gene regions. The multi-gene region loci weighted method is the most complex but also the most effective. Its power simulation results are much higher than the unweighted single gene region analysis method and the false positive ratio is much lower than the single gene region analysis method. For the SLoS method, the simulation results of the common multi-gene region method are only slightly better than the common single-gene region analysis method, which also shows that even the simplest multi-gene region analysis method can effectively improve the power of the analysis compared with the single-gene region analysis method. In general, the simulation results of the Step method and LW-Step method are better methods for the associated analysis of multi-gene region. By modified freedom degree of test statistic F, the multi-SLoS, W-SLoS and LW-SLoS are feasible for the associated analysis of multi-gene region. Compared with the rare variant multi-gene region, the associated analysis result of the common variant multi-gene region is better than using the multi-gene region analysis method.
The multi-gene region loci weighted method is a further expansion and extension of the weighting idea of Belonogova et al. [32]. However, there are some differences between our work and Belonogova's paper: firstly, the weighted idea was not applied to a multigene region in his paper; secondly, the coefficient of functional data is estimated by the smooth method in our paper and the least square method in Belonogova's paper.
The Fourier basis function is selected to fit the genotype data when the functional expansion is carried out. The reason as to why we chose the Fourier basis function is that some studies have compared the Fourier basis to the spline basis, achieving similar results (the previously cited papers on functional gene-association analysis have compared it in their papers). However, some people question how genes can be represented by periodic Fourier bases. Perhaps they both extracted the same amount of information for the gene regions they wanted to summarize in their way, and after a lot of people compared the results of the two bases, there was only a very small difference. Even though the Fourier basis function might look better in some cases, it is not good enough for the authors of those papers to conclude that the Fourier basis function is a better choice. Usually, both of the basis functions are used, so one can choose one of two types of basis functions. The reason why the Fourier basis is selected in this paper is that the Fourier basis only needs to determine the number of basis functions, while the spline basis not only needs to select the number of basis functions but also the order of the basis function, so the selection of the Fourier basis function is much simpler. It must choose the Fourier basis for the Step method and the LW-Step method.
Regarding the selection of compression parameters and smoothing parameters, the paper does not necessarily choose the weights that can give the method the best performance, but it basically selects the parameters according to the most common standards and methods (such as 10-fold cross-validation, etc.). The specific parameter selection strategy can be: Firstly, determine the selection criteria of parameters, then determine the approximate range of compression parameters and smoothing parameters according to the method and the actual situation, finally, the computer program is used to screen out the optimal parameters. This process must be repeated in the processing of actual data, because the genetic region composition of the actual data is constantly changing, and so it requires that the parameters change accordingly. However, in the simulation, the same situations in the genetic regions are made up the same way. Therefore, in order to save computing resources, a suitable parameter is directly selected in this paper.
One might ask: why not just analyze it as a larger gene region? Instead, we analyze it in the way of this multi-gene region, and the functional data can do this only by increasing the number of nodes and the number of basis functions. One reason is that the large gene regions can be divided into smaller gene regions, and then the multi-gene region method can be used to get more detailed test results. Another reason is that, as mentioned in the previous article, due to the consideration of the interrelationship between different gene regions, the multi-gene region method can better identify some regions with relatively weak effects and have higher power. So, compared to single-gene region testing, multigene region testing can detect gene regions where the effect is smaller or where the effect overlaps with that of other regions.
We focused on independent SNPs for common variants. Rare variants come from the rare haploid dataset (the R language SKAT package [33] contains a set of such data produced by simulating allele frequency and linkage imbalance information in European populations) with a length of 200 kb. However, we also simulate linkage disequilibrium in common variant multi-gene region of Scenario I based on the common multi-gene region method and multi-gene region weighted method. When the r 2 measure of linkage disequilibrium is between 0.25 and 0.64, the power of the linkage disequilibrium based on the SLoS and Multi-SloS methods is higher than that of linkage equilibrium, but the false positive rates also increase significantly. The power of the linkage disequilibrium based on the Step method decreases slightly, and the false positive rates remain largely unchanged. The power and false positive rates of the linkage disequilibrium based on the FLM method does not change significantly. This indicates that the linkage disequilibrium among the gene loci causes SLoS and Multi-SLoS methods to more easily misidentify nonassociated gene regions. The association analysis between gene region and quantitative trait is susceptible to linkage disequilibrium of gene loci for SLoS and Multi-SLoS methods. Then, the association analysis of gene regions for quantitative traits was unstable for SLoS and Multi-SLoS methods, while Step and FLM methods are more stable. The reason may be the parameter estimation of SLoS and Multi-SLoS method. Furthermore, we compare the simulation results of forward selection and backward selection and find that the power of the forward selection method is slightly higher than that of the back selection method, but the false positive rates of the forward selection method is far greater than that of the back selection method. Finally, we study the distribution of allele frequencies of gene loci. The results show that the power of following a normal distribution is higher than that of following a uniform distribution, which may be due to the fact that the MAF of the normal distribution is mostly smaller than that of the uniform distribution. According to the effect function expression, the effect size of the normal distribution is larger at this time and can make the locus easier to detect.
The model we considered is an ideal state, as it is only a study of basic assumptions, without considering group structure and population with relatives, and there is no missing genotype. In practice, it is inevitable that there will be missing genotypes. At this time, we can use statistical methods to fill in the missing data, and then convert the discrete genetic data into continuous functions. In future research, we will consider incorporating models of group structure and population with relatives.