Analytical method, system and storage medium capable of controlling kinship correlation in samples

Abstract

The present invention discloses an analysis method, system and storage medium capable of controlling kinship correlation in a sample, comprising the following steps: S1: acquiring sample data; S2: constructing a generalized linear model based on the sample data, estimating the parameters of a fitted constraint model and calculating residuals; S3: dividing the residuals of the fitted constraint model into outliers and non-outliers, calculating and storing the homology sharing probability between any two members in a family where the number of members is greater than 1 and the sample data corresponding to the outlier is located; S4: dividing all sites in the whole genome into a certain number of intervals according to the minor allele frequency, calculating and storing the joint probability of the genotype vectors of each family member at the endpoints of each interval; S5: calculating the statistical p-value for each site to be tested using a hybrid test strategy combining a normal distribution approximation method with an empirical saddle point approximation method; the analysis method avoids a large number of false positive or false negative results.

Description

Technical Field

[0001] The present invention relates to the field of biological technology, and in particular to an analysis method, system and storage medium capable of controlling kinship correlation in a sample. Background Art

[0002] Genome-wide association studies (GWAS) often involve screening the entire human genome for sequence variants, such as single-nucleotide polymorphisms (SNPs), that are significantly associated with specific traits. In recent years, with the advancement and application of gene sequencing and electronic health record technologies, biobanks have accumulated vast amounts of diverse health data, becoming a crucial source of data for GWAS studies of complex traits. For example, the UK Biobank has collected health data from 500,000 British individuals, making them available for analysis by researchers worldwide. This data includes multidimensional, cross-scale data such as gene sequencing, magnetic resonance imaging, clinical indicators, and lifestyle data. Complex traits can provide more information than simpler traits. For example, longitudinal data can reflect changes in an individual's health status over time, providing richer information than conventional cross-sectional data. Using rapid and accurate methods to extract useful information from complex traits is crucial for discovering causal genes, exploring the mechanisms of complex diseases, and implementing precision medicine and personalized prevention.

[0003] Large biobanks contain large amounts of data, numerous confounding factors, and complex data formats, necessitating rapid, accurate, and universal methods for analysis. GWAS are often applied to hundreds of thousands of samples and millions of SNPs, placing a significant computational burden on association analysis methods. Common testing methods, such as the likelihood ratio test and the Wald test, require fitting an unconstrained model under the alternative hypothesis for each locus to be tested. This requires millions of fits across the entire genome, resulting in a very long computational time. The score test, on the other hand, only requires fitting a constrained model under the null hypothesis once across the entire genome, greatly improving computational efficiency. However, unbalanced phenotypic and genotypic distributions can violate the normal distribution assumption of the score statistic, resulting in an excessively high false positive rate in the statistical test. The saddlepoint approximation (SPA) method utilizes information from higher-order moments of the score statistic to obtain a more accurate distribution approximation, effectively controlling the Type I error rate.

[0004] Kinship is an important confounding factor in large biobanks. For example, about one-fifth of the samples in the UK Biobank are kin-related. If the kinship problem cannot be well controlled, a large number of false-positive sites will be generated. Eliminating data with kinship in large biobanks will significantly reduce the sample size and statistical power. Linear mixed models are often used to control kinship problems. The basic idea is to introduce a random effect term into the model and assume that the random effect term obeys a multivariate Gaussian distribution with a mean of zero vector and a covariance matrix of the genetic relationship matrix (GRM). However, it is very difficult to apply the saddle point approximation method and linear mixed models to the analysis of complex traits.

[0005] Complex traits in large biobanks manifest in a variety of forms, including quantitative, dichotomous, time-series, multi-category, and longitudinal data. In recent years, numerous GWAS analysis methods have been developed for analyzing quantitative, dichotomous, time-series, and multi-category traits, such as fastGWA, SAIGE, GATE, and POLMM. These methods all perform genome-wide association analysis based on the score test, apply linear mixed models to control for kinship within the sample, and employ saddle point approximations to accurately determine the distribution of the score statistic, enabling rapid and accurate calculation of statistical p-values. However, suitable GWAS analysis methods for longitudinal data and more complex data types are still lacking.

[0006] The TrajGWAS method is a recently proposed GWAS method suitable for analyzing longitudinal trait data from large biobanks. In addition to using between-subject variability (BS variability) to model the mean level of individual traits, it also uses within-subject variability (WS variability) to capture the magnitude of temporal fluctuations in individual traits, a key risk factor for many diseases. Furthermore, the TrajGWAS method employs an empirical saddlepoint approximation to accurately obtain the distribution of the score statistic, enabling it to control the Type I error rate when testing associations between both within-subject and between-subject variability in longitudinal traits and single-nucleotide polymorphisms. However, the empirical saddlepoint approximation is only applicable to the analysis of independent samples and will generate a large number of false-positive results in the analysis of kin-related samples. Therefore, the TrajGWAS method cannot analyze data from kin-related samples, resulting in a high number of false-positive results. Excluding these related samples significantly reduces statistical power.

[0007] In summary, there is an urgent need for an analytical method that is suitable for complex traits and can control the kinship correlation in samples to avoid a large number of false positive or false negative results. Summary of the Invention

[0008] Based on the technical problems existing in the background technology, the present invention proposes an analysis method, system and storage medium that can control the kinship correlation in samples, which effectively controls the first type error rate and avoids a large number of false positive or false negative results; it can be applied to analyze samples with kinship correlation to obtain higher statistical power.

[0009] The present invention proposes an analytical method for controlling kinship correlation in samples, comprising the following steps:

[0010] S1: Obtain sample data, wherein the sample data includes phenotypic data, genotypic data, and confounding factor data;

[0011] S2: constructing a generalized linear model based on the sample data, ignoring the random effect terms in the generalized linear model, fitting the generalized linear model under the assumption that the genetic effect is zero to obtain a fitted constrained model, estimating the parameters of the fitted constrained model and calculating the residuals;

[0012] S3: Divide the residuals of the fitted constraint model into outliers and non-outliers, and divide the sample data into K preFor each family, the number of family members greater than c and containing sample data corresponding to the outlier is divided into K families containing sample data corresponding to the outlier. The probability of homology sharing between two members in the family with more than 1 and containing sample data corresponding to the outlier is calculated and stored.

[0013] S4: Divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency, and calculate and store the joint probability of the genotype vectors of each family member at the endpoints of each interval based on the homologous sharing probability and the Chow-Liu algorithm;

[0014] S5: For each site to be tested, estimate the minor allele frequency and calculate the score statistic. Use a hybrid test strategy that combines the normal distribution approximation method with the empirical saddle point approximation method to calculate the statistical p-value to achieve genome-wide association analysis.

[0015] Furthermore, in step S2, it specifically includes:

[0016] Assume that the number of sample data is n, let X represent an n×p dimensional confounding factor matrix, G represent the genotype vector, and Y represent a phenotypic data;

[0017] Calculate the linear prediction term η=Xβ based on X, G, and Y X +Gβ G +b, treat the linear prediction term as a generalized linear model, where β X represents the p-dimensional coefficient vector of the confounding factor X, β G represents the genetic effect of genotype G, b represents the random effect term used to describe the relationship between individuals, b ~ N (0, τ∑), where ∑ is the n × n dimensional genetic relatedness matrix (GRM), and τ is the additive genetic variance;

[0018] Ignoring the random effect term b, the genetic effect β G Under the assumption that is zero, a generalized linear model is fitted to obtain a fitted constrained model, and the fitted constrained model is used to obtain the estimated values ​​of the parameters under the constrained conditions and to calculate the residuals of the fitted constrained model;

[0019] Based on the estimated value, we can calculate the G The statistic S is zero.

[0020] Furthermore, in step S3, it specifically includes:

[0021] Preselect a positive integer γ and let α 0. and α 0.They represent the 25% and 75% quantiles of the residual R obtained by fitting the constrained model, and the interquartile range is defined as IQR = α 0.75 -α 0.25 ;

[0022] R i >α 0.75 +IRQ·γ or R i <α 0.25 The residual corresponding to -IRQ·γ is regarded as an outlier, and the other residuals are regarded as non-outliers;

[0023] The samples are divided into K groups based on the affinity coefficient of the sample data. pre A family, so that the kinship coefficient between family members is not 0, and the kinship coefficient between different family members is equal to 0;

[0024] Preselect a positive integer c, split the families with more than c members and containing sample data corresponding to outliers, and obtain K families containing sample data corresponding to outliers;

[0025] Based on any K families whose number of members is greater than 1 and contains the sample data corresponding to the outlier, set the i-th family to contain members 1,...,K i ,1≤i≤K, the corresponding residual is make represents the homologous sharing probability that the jth and kth members in the i-th family share 0, 1, and 2 genetic variants at the same locus, respectively;

[0026] Under the assumption of Hardy-Weinberg equilibrium, the genotype variables of each member follow a binomial distribution;

[0027] Using the information of s sites in the whole genome, we can obtain the kinship coefficient ρ for the jth and kth members. jk The probability of sharing 0 genetic variation at the same site with members j and k is λ jk Moment Estimates

[0028] Moment-based estimation Calculate the probability of homology sharing Moment Estimates

[0029] The moment estimate and moment estimates The calculation formula is as follows:

[0030]

[0031]

[0032]

[0033]

[0034] in, represents the estimated value of the kinship coefficient between members j and k, represents the estimated probability of homologous sharing between members j and k sharing 0 genetic variants at the same site, G jl and G kl denote the genotype of the lth site of the jth and jth members, respectively. represents the weight of the genotype at the lth site, μ l represents the minor allele frequency of the lth site, represents μ l Estimates.

[0035] Furthermore, a positive integer c is pre-selected, and families with a number of members greater than c and containing sample data corresponding to outliers are segmented, specifically including:

[0036] Preselect a positive integer c;

[0037] Assume that any family i with more than c members and containing outliers contains members The corresponding residual is

[0038] Calculate the interaction parameter |ρ between any two members jk R j R k |, where ρ jk represents the kinship coefficient between the jth and kth members, R j and R k represents the corresponding residual;

[0039] The interaction parameter |ρ jk R j R k Arrange them in ascending order, and delete the related member pairs in sequence until the number of family members after the split is no more than c, and obtain K families that contain sample data corresponding to the outliers.

[0040] Furthermore, in step S4, it specifically includes:

[0041] Divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency;

[0042] Based on the fact that any family member number is greater than 2 and contains the sample data corresponding to the outlier, it is assumed that any family i contains members 1,...,K i ,K i>2, let G i Represents the genotype vector corresponding to the family members Let Pr(G i ) represents the joint probability of the genotype vector of the family member;

[0043] Based on the genotype variable G of the jth and kth members j and G k Mutual information I(G j ,G k ), minimize the joint probability distribution Pr(G i ) and the approximate joint probability distribution The relative entropy of Minimizing the KL divergence is equivalent to calculating the value of I(G j ,G k ),1≤j <k≤K i Construct (K i -1) Maximum spanning tree of edges;

[0044] Using Prim's algorithm according to (K i (K i -1) / 2) mutual information I(G j ,G k ),1≤j <k≤K i , construct (K i -1) Maximum spanning tree of edges;

[0045] Based on the Chow-Liu algorithm, according to (K i -1) The maximum spanning tree of edges obtains the optimal approximate joint probability distribution Approximate joint probability distribution As the genotype vector G of the family member i The joint probability of

[0046] Repeat the above operation at each interval endpoint to calculate the joint probability of the genotype vectors of other family members, and store all the joint probabilities to obtain the joint probability of the genotype vectors of each family member at each interval endpoint corresponding to all sites in the whole genome.

[0047] Furthermore, the I(G j ,G k ) is calculated as follows:

[0048]

[0049]

[0050]

[0051] Among them, m l Indicates from 1 to K i Rearrangement, 1≤l≤K i , m h(l) represents the rearrangement from 0 to l-1, 0≤h(l) <l, Indicates that at a given mth h(l) Member's genotype Under the condition m l Member's genotype The conditional probability of I(G j ,G k ) represents the genotype variable G of the jth and kth members j and G k The mutual information of H(G l ) is the genotype variable G of the lth member l The information entropy of H(G i ) represents the genotype vector of the family member The joint entropy of g j ,g k ∈{0,1,2} represents the genotype variable G of the jth and kth members respectively j and G k The possible values ​​of P(g j ),P(g k ) represent the genotype G of the jth and kth members respectively j and G k The observed value is g j ,g k The probability, P(g j ,g k ) represents the genotype variable G of the jth and kth members j and G k The observed value is g j ,g k The joint probability of .

[0052] Furthermore, in step S5, it specifically includes:

[0053] Assume that the sample genotype vector of a certain site to be tested is G=(G1,G2,…,G n ) T , the corresponding residual is R=(R1,R2,…,R n ) T ;

[0054] Under the assumption of Hardy-Weinberg equilibrium, the genotype variables of each member follow a binomial distribution, and the genetic effect β is calculated. G The score statistic S under the assumption that it is zero;

[0055] (a) When using the normal distribution approximation method, the details are as follows:

[0056] Calculate the expectation and variance of the score statistic S, and calculate the statistical p-value based on the expectation and variance:

[0057]

[0058]

[0059]

[0060]

[0061] in, is the observed value of S, Φ(·) is the cumulative distribution function of the standard normal distribution, Var(G i ) represents the genotype variable G of the i-th individual i The variance of Cov(G i ,G j ) represents the genotype variable G of the i-th and j-th individuals i and G j The covariance of g Represents the genotype variable G i The standard deviation of μ represents the minor allele frequency of the tested locus, and ρ ij represents the kinship coefficient between members i and j, R i and R j Denotes the corresponding residual; let As an estimate of μ, As σ g Estimates, As an estimate of Var(S);

[0062] (b) Use the saddle point approximation method as follows:

[0063] Decompose the score statistic S according to the kinship coefficient of the sample data and the corresponding residual to obtain the score statistic S of the family corresponding to the sample data of the outlier residual. outliers and the score statistic S of the family corresponding to the sample data without outlier residuals non-outliers ;

[0064] According to S outLiers and S non-outliers Calculate the moment generating function M of the score statistic S (t);

[0065] The moment generating function M S (t) After accumulation, the moment function K is generated S (t);

[0066] Based on the moment function K S (t) Calculate the statistical p value;

[0067] in:

[0068] S=S outliers +S non-outliers

[0069]

[0070]

[0071]

[0072]

[0073]

[0074]

[0075]

[0076]

[0077]

[0078] Among them, S i represents the score statistic of the family containing the sample data corresponding to the outlier, represents the genotype of the jth member of the i-th family containing the sample data corresponding to the outlier, Indicates S non-outliers The moment generating function of Indicates S outliers The moment generating function, M Si (t) represents the moment generating function of the score statistic of the family containing the sample data corresponding to the outlier, E(S non-outliers ) indicates S non-outliers Expectations, R non-outliers represents the residual corresponding to the members of the family excluding the sample data corresponding to the outlier, Represents R non-outliers The transpose of μ represents the minor allele frequency, Var(S non-outliers ) indicates S non-outliers The variance of represents the variance of the genotype variable at the site to be tested, ∑ non-outliers Represents the genetic correlation matrix between members of a family excluding sample data corresponding to outliers, where the jth and kth elements represent the kinship coefficient ρ of the jth and kth members jk, K′ S (t) and K″ S (t) represents the moment function K S The first and second order derivatives of (t), the saddle point ζ satisfies in, represents the observed value of S. According to the Barndorff-Nielsen saddle point approximation formula, in Φ(·) represents the cumulative distribution function of the standard normal distribution, and then the two-sided statistical p value is obtained in yes Estimates, yes Estimates.

[0079] Further, The calculations vary slightly depending on the number of family members, including:

[0080] If the family i containing the sample data corresponding to the outlier has only one member i1, and the genotype variable of member i1 follows the binomial distribution, we get Among them, R i1 represents the residual corresponding to member i1, μ represents the minor allele frequency of the site to be tested;

[0081] If the family i containing the sample data corresponding to the outlier contains two members i1 and i2, the genotype variables of members i1 and i2 follow the binomial distribution, and the homologous sharing probability estimate of members i1 and i2 sharing 0, 1, and 2 genetic variations at the tested site is and the corresponding S i Moment generating function

[0082]

[0083]

[0084]

[0085] The sum of the products gives in Respectively represent the residuals corresponding to members i1, i2, S represents the case where members i1 and i2 share 0, 1, and 2 genetic variations at the same site, respectively. i The moment generating function of , μ represents the minor allele frequency of the site to be tested;

[0086] If the family with the outlier corresponding to the sample data of the i-th family contains more than two members According to the estimated value of the joint probability of the genotype vectors of family members at the endpoints of each minor allele frequency interval, when the minor allele frequency μ of the site to be tested falls into a certain minor allele frequency interval, the joint probability of the genotype vectors of the family members is It is approximated by the linear weighted average of the estimated values ​​of the joint probability at the two end points of the interval, S i Moment generating function

[0087]

[0088] An analysis system capable of controlling kinship correlation in a sample comprises an acquisition module, a construction calculation module, a homology sharing probability calculation module, a joint probability calculation module and a statistical value calculation module;

[0089] The acquisition module is used to acquire sample data, and the sample data includes phenotypic data, genotypic data and confounding factor data;

[0090] The construction calculation module is used to construct a generalized linear model based on the sample data, ignore the random effect terms in the generalized linear model, fit the generalized linear model under the assumption that the genetic effect is zero to obtain a fitted constraint model, estimate the parameters of the fitted constraint model and calculate the residual;

[0091] The homology sharing probability calculation module is used to divide the residuals of the fitted constraint model into outliers and non-outliers, and divide the samples into K pre For each family, the number of members is greater than c and the families containing sample data corresponding to the outlier are divided to obtain K families containing sample data corresponding to the outlier. The probability of homology sharing between two members in the family with more than 1 and containing sample data corresponding to the outlier is calculated and stored.

[0092] The joint probability calculation module is used to divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency, and calculate and store the joint probability of the genotype vectors of each family member at the endpoints of each interval based on the homologous sharing probability and the Chow-Liu algorithm;

[0093] The statistical value calculation module is used to estimate the minor allele frequency and calculate the score statistic for each site to be tested, and calculate the statistical p-value using a hybrid test strategy that combines the normal distribution approximation method with the empirical saddle point approximation method to achieve genome-wide association analysis.

[0094] A computer-readable storage medium having a plurality of programs stored thereon, wherein the plurality of programs are used to be called by a processor and execute the method for analyzing kinship correlation in controllable samples as claimed in any one of claims 1 to 8.

[0095] Those skilled in the art will understand that all or part of the steps of implementing the above-mentioned method embodiment can be completed by hardware related to program instructions, and the aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it executes the steps of the above-mentioned method embodiment; and the aforementioned storage medium includes: ROM, RAM, disk or optical disk, etc. Various media that can store program codes.

[0096] The advantages of the analysis method, system and storage medium for controlling kinship correlation in samples provided by the present invention are: the analysis method, system and storage medium for controlling kinship correlation in samples provided in the structure of the present invention can well control the first type error rate; it can be applied to analyze samples with kinship correlation to obtain higher statistical power, specifically: fitting a fitting constraint model that does not include genotypes, that is, fitting a constraint model and model parameters under the assumption that the genetic effect is zero and calculating the residuals; dividing the residuals of the fitting constraint model into outliers and non-outliers, and dividing the samples into K groups based on the kinship coefficient of the sample data. pre The invention provides a method for analyzing the genetic association of a whole genome in a fast and accurate manner. The invention provides a method for analyzing the genetic association of a whole genome in a fast and accurate manner. The method comprises the following steps: first, dividing the families with a number of members greater than c and containing sample data corresponding to the outliers into K families containing sample data corresponding to the outliers; second, calculating and storing the homology sharing probability between two members in the families with a number of members greater than 1 and containing sample data corresponding to the outliers; third, dividing all sites in the whole genome into a certain number of intervals according to the minor allele frequency (MAF); fourth, calculating and storing the joint probability distribution of the genotype vectors of each family at the endpoints of each interval using the Chow-Liu algorithm; fifth, estimating the MAF of each site to be tested and calculating the score statistic; fourth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; fifth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; fifth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; fifth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; sixth ... BRIEF DESCRIPTION OF THE DRAWINGS

[0097] Figure 1 It is a schematic diagram of the process of the present invention;

[0098] Figure 2 Schematic diagram of a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method;

[0099] Figure 3This is the genetic pedigree diagram used to simulate samples with kinship in numerical simulation experiments. The families in population A have genetic relationships as shown on the left, and the families in population B have genetic relationships as shown on the right.

[0100] Figure 4 This is a numerical simulation result of the first type error rate in the longitudinal data phenotyping analysis for the case where the genetic effect of the site on both intra-individual variation and inter-individual variation is 0, where A represents the test β G =0, B represents the test τ G =0, the horizontal axis represents the type of variant site (common variant site and rare variant site), and the vertical axis represents the first type error rate at the corresponding significance level; considering two significance levels α = 5×10 -7 and α=5×10 -5 Consider three populations: population A: small family structure population, population B: large family structure population and population C: unrelated population; consider comparing three methods: a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method (SPA GRM ), Normal distribution approximation method (Norm GRM ) and TrajGWAS methods;

[0101] Figure 5 This is a numerical simulation result of statistical power in longitudinal data phenotyping, where A represents the test β G =0, B represents the test τ G =0, the horizontal axis represents the hybrid test strategy (SPA) combining the two methods, namely the TrajGWAS method and the normal distribution approximation method and the saddle point approximation method. GRM ), the vertical axis represents the chi-square value corresponding to the p-value; two types of variant sites are considered, namely common variant sites and rare variant sites; three populations are considered, namely population A: small family structure population, B: large family structure population and C: unrelated population. DETAILED DESCRIPTION

[0102] The technical solutions of the present invention are described in detail below through specific embodiments. Numerous specific details are set forth in the following description to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0103] like Figures 1 to 5 As shown, the analysis method for controlling kinship correlation in samples proposed by the present invention includes the following steps S1 to S5.

[0104] S1: Obtain sample data, wherein the sample data includes phenotypic data, genotypic data, and confounding factor data;

[0105] The confounding factor data is a series of data including confounding factors such as age and gender. Assuming that the phenotypic data is longitudinal data, in the above step S1, it is assumed that the number of individual samples participating in the study is n, m i represents the number of longitudinal data observations for the i-th individual, Represents the total number of longitudinal observations for all samples. X, Z, and W represent a m×p, m×q, and m×l dimensional confounding factor matrix (including confounding factors such as age and gender), respectively. X, Z, and W can contain the same confounding factors, and G represents the genotype vector (G i =0,1,2), Y represents the observed value of the longitudinal phenotype.

[0106] S2: constructing a generalized linear model based on the sample data, ignoring the random effect terms in the generalized linear model, fitting the generalized linear model under the assumption that the genetic effect is zero to obtain a fitted constraint model, estimating the parameters of the fitted constraint model and calculating the residual, specifically including steps S21 to S23:

[0107] S21: Assume that the number of sample data is n, let X represent an n×p-dimensional confounding factor matrix, G represent the genotype vector, and Y represent a certain phenotypic data;

[0108] The X in step S21 has the same meaning as that in step S1, except that it is used to represent an n×p-dimensional confounding factor matrix in step S12, while X in S1 is used to represent an m×p-dimensional confounding factor moment. The specific values ​​of n, m, and p are selected according to actual conditions, and thus the two Xs are essentially the same.

[0109] S22: Calculate the linear prediction term η=Xβ based on X, G, and Y X +Gβ G +b, treat the linear prediction term as a generalized linear model, where β X represents the p-dimensional coefficient vector of the confounding factor X, β G represents the genetic effect of genotype G, b represents the random effect term used to describe the relationship between individuals, and it is assumed that b~N(0,τ∑), where ∑ is the n×n-dimensional genetic relationship matrix (GRM) and τ is the additive genetic variance.

[0110] The above linear prediction term is an overall summary formula expression. The b in the linear prediction term is also a general term. The following will explain it in detail:

[0111] For the longitudinal data phenotypes, the following mixed-effects multiple location scale model was fit as a generalized linear model:

[0112]

[0113]

[0114] Among them, y ij represents the j-th longitudinal phenotypic observation value of the i-th individual; x ij ,z ij ,w ij Represents a p×1, q×1, l×1 dimensional confounding factor vector (including confounding factors such as age and gender), where x ij ,z ij ,w ij May contain the same confounding factors, G i represents the genotype of the i-th individual; β, τ, respectively represent x ij ,w ij The corresponding coefficient, β G ,τ G Indicates G i The corresponding coefficient of G ,τ G Respectively represent the genetic effects of the tested loci on between-subject variability and within-subject variability, ε ij represents the error term of the j-th longitudinal phenotypic observation of the i-th individual, Represents x ij ,z ij ,w ij The transpose of ε ij ,ω i ,b i , represent random effect terms respectively.

[0115] Assuming that the random effect term (γ i T ,ω i ) T Obey a mean vector of 0 q+1 , the covariance matrix is Multivariate normal distribution; assuming the random effect term b i , Each obeys a mean vector of 0 n , the covariance matrix is ​​σ 2∑ and The multivariate normal distribution of ∑ is the n×n dimensional genetic relationship matrix (GRM), σ 2 and denote the additional genetic variance of the corresponding equations, and assume that the error term ε in equation (1) ij The variance of satisfies equation (2).

[0116] S22: Ignore the random effect term b, and in the genetic effect β G Under the assumption that is zero, a generalized linear model is fitted to obtain a fitted constrained model, and the fitted constrained model is used to obtain the estimated values ​​of the parameters under the constrained conditions and to calculate the residuals of the fitted constrained model;

[0117] In order to test the significance of the genetic effect of gene-longitudinal traits, the above mixed-effect multi-location scale model was fitted and calculated, specifically: the random effect term b in formulas (1) and (2) i and Using the WiSER parameter estimation method, under the assumption that the genetic effect is zero, H0:β G =τ G = 0, and the mixed-effect multi-location scale model is fitted and the parameters are calculated. The WiSER parameter estimates are:

[0118]

[0119]

[0120]

[0121]

[0122]

[0123] is the covariance matrix Var(Y i ) and use τ,∑ γ The least squares estimator of Each subsequent iteration and Until Converge to obtain the WiSER parameter estimates for the fitted constrained model ∑ is the n×n dimensional genetic relationship matrix (GRM), and τ is the additive genetic variance.

[0124] S23: Calculate the value based on the estimated value in βG The statistic s when it is zero;

[0125] Fit a constrained model, estimate parameters, and calculate the test beta under the assumption that the genetic effect is zero G =0 and τ G =0, the score statistics are:

[0126]

[0127]

[0128] in Can be obtained from the fitted constraint model, defined as the fitted constraint model and Test β G =0 and τ G =0 residual.

[0129] S3: Divide the residuals of the fitted constraint model into outliers and non-outliers, and divide the sample data into K pre families, splitting the families with more than c family members and containing sample data corresponding to the outlier to obtain K families containing sample data corresponding to the outlier, and calculating and storing the homology sharing probability between any two members of the families with more than 1 family member and containing sample data corresponding to the outlier, specifically including steps S31 to S38:

[0130] S31: Preselect a positive integer γ, for example, γ = 1.5, and let α 0.25 and α 0.75 They represent the 25% and 75% quantiles of the residual R obtained by fitting the constrained model, and the interquartile range (IQR) is defined as IQR = α 0.75 -α 0.25 ;

[0131] definition is the test β obtained from fitting the constrained model G =0(τ G =0), then test β G =0(τ G =0) is recorded as S = R T G.

[0132] S32: R i >α 0.75 +IRQ·γ or R i <α 0.25 The residuals corresponding to -IRQ·γ are regarded as outliers, and the other residuals are regarded as non-outliers.

[0133] S33: Divide the sample into K groups based on the affinity coefficient of the sample data pre A family, so that the kinship coefficient between family members is not 0, and the kinship coefficient between different family members is equal to 0;

[0134] The (j,k)th element in the genetic correlation matrix GRM is the kinship coefficient ρ of individuals j and k jk .

[0135] S34: Preselect a positive integer c, such as c=5, and segment the families whose number of members is greater than c and contain sample data corresponding to the outliers, to obtain K families containing sample data corresponding to the outliers;

[0136] To improve computational efficiency, consider segmenting families whose number of members is greater than c and that contain samples corresponding to outliers. To maintain computational accuracy, use a greedy algorithm to segment families that meet the conditions into multiple families whose number of members is no more than c. The specific approach is as follows: Steps S341 to S343:

[0137] S341: Assume that any family i with a number of members greater than c and containing an outlier corresponding to the sample data contains members 1≤i≤K pre , the corresponding residual is

[0138] S342: Calculate the interaction parameter |ρ between any two members jk R j R k |, where ρ jk represents the kinship coefficient between the jth and kth members, R j and R k represents the corresponding residual;

[0139] S343: The interaction parameter |ρ jk R j R k Arrange them in ascending order, and delete the related member pairs in sequence until the number of family members after the split is no more than c, and obtain K families that contain sample data corresponding to the outliers.

[0140] Through steps S341 to S343, not only the analysis and calculation efficiency is improved, but also the calculation accuracy is maintained.

[0141] S35: Based on any K families whose number of members is greater than 1 and contains sample data corresponding to the outlier, set the i-th family to contain members 1, ..., K i ,1≤i≤K, the corresponding residual is make represents the homologous sharing probability that the jth and kth members in the i-th family share 0, 1, and 2 genetic variants at the same locus, respectively;

[0142] S36: Under the assumption of Hardy Weinberg equilibrium, the genotype variables of each member follow a binomial distribution, that is, G i ~B(2,μ), where μ represents the minor allele frequency, precomputed:

[0143]

[0144]

[0145] where ρ jk represents the kinship coefficient between members j and k, λ jk represents the homologous sharing probability that the jth and kth members share 0 genetic variants at the same site, and μ represents the minor allele frequency (MAF) of this site.

[0146] S37: Use the information of s sites in the whole genome to obtain the kinship coefficient ρ for the jth and kth members jk The probability of sharing 0 genetic variation at the same site with members j and k is λ jk Moment Estimates

[0147] In practice, the information of s sites is used to obtain the jk and λ jk Moment estimate of :

[0148]

[0149]

[0150]

[0151] in, represents the estimated value of the kinship coefficient between members j and k, represents the estimated probability of homologous sharing between members j and k sharing 0 genetic variants at the same site, G jl and G kl denote the genotype of the lth site of the jth and kth members, respectively. represents the weight of the genotype at the lth site, μ l represents the minor allele frequency of the lth site, represents μ l Estimates.

[0152] S38: Moment-based estimation Calculate the probability of homology sharing Moment Estimates

[0153] The following linear equations are obtained by solving Moment estimate of :

[0154]

[0155] For example, the probability of a father and son (father and daughter / mother and son / mother and daughter) sharing (0, 1, 2) genetic materials is p. (0) =0,p (1) =1,p (2) =0, kinship coefficient For example, the probability of a pair of brothers (sisters / brothers / sisters) sharing 0, 1, and 2 genetic materials is Kinship coefficient

[0156] S4: Divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency, and calculate and store the joint probability of the genotype vectors of each family member at the endpoints of each interval based on the homologous sharing probability and the Chow-Liu algorithm, specifically including steps S41 to S46:

[0157] S41: Divide all loci in the whole genome into a certain number of intervals according to the minor allele frequency (MAF);

[0158] For example, MAF∈(0.0001,0.5], then this range is divided into a certain number of intervals, such as 10 -4 ,5×10 -4 ,10 -3 ,5×10 -3 ,0.01,0.05,0.1,0.2,0.3,0.4,0.5.

[0159] S42: Based on any family with more than 2 members and containing the sample data corresponding to the outlier, set any family i to contain members 1, ..., K i ,K i >2, let G i Represents the genotype vector corresponding to the family members Let Pr(G i ) represents the joint probability of the genotype vector of the family member;

[0160] The core idea of ​​the Chow-Liu algorithm is to use (n-1) second-order conditional distributions to optimally approximate the joint probability distribution of n random variables.

[0161] Pr(G i) is approximately:

[0162]

[0163] Among them, m l Indicates from 1 to K i Rearrangement, 1≤l≤K i , m h(l) represents the rearrangement from 0 to l-1, Indicates that at a given mth h(l) Member's genotype Under the condition m l Member's genotype The conditional probability of

[0164] S43: Genotype variable G based on the jth and lth members j and G l Mutual information I(G j ,G l ), minimize the joint probability distribution Pr(G i ) and the approximate joint probability distribution The relative entropy of Minimizing the KL divergence is equivalent to I(G j ,G l ) Construct k i Maximum spanning tree of edges;

[0165] To obtain Middle m j and the corresponding m h(j) , minimize Pr(G i )and The relative entropy, or KL divergence (Kullback-Leibler divergence):

[0166]

[0167]

[0168] Among them, I(G j ,G k ) represents the genotype variable G of the jth and kth members j and G k The mutual information of H(G l ) is the genotype variable G of the lth member l The information entropy of H(G i ) represents the genotype vector of the family member The joint entropy of g j ,g k∈{0,1,2} represents the genotype variable G of the jth and kth members respectively j and G k The possible values ​​of P(g j ),P(g k ) represent the genotype variables G of the jth and kth members respectively j and G k The observed value is g j ,g k The probability, P(g j ,g k ) represents the genotype variable G of the jth and kth members j and H k The observed value is h j ,g k The joint probability of .

[0169] Since KL divergence does not depend on ∑H(G l ) and H(G i ), only depends on ∑I(G j ,G k ), so minimizing the KL divergence is equivalent to according to I(G j ,G k ) to construct a maximum spanning tree.

[0170] S44: Using Prim's algorithm according to (K i (K i -1) / 2) mutual information I(G j ,G k ),1≤j <k≤K i , construct (K i -1) Maximum spanning tree of edges;

[0171] S45: Based on the Chow-Liu algorithm, determine (K i -1) edges, thereby determining m l and the corresponding m h(l) , and obtain the optimal joint probability approximation Approximate joint probability distribution as the joint probability of the genotype vectors of the family members.

[0172] S46: Due to The calculation of depends on the minor allele frequency (MAF) of the locus to be analyzed. Therefore, the above operation is repeated at each endpoint of the interval to calculate the joint probability of the genotype vectors of other families, and all the joint probabilities are stored to obtain the joint probability of the genotype vectors of each family member at each endpoint of each interval corresponding to all sites in the whole genome.

[0173] S5: For each site to be tested, estimate the minor allele frequency and calculate the score statistic, and use a hybrid test strategy that combines the normal distribution approximation method with the empirical saddle point approximation method to calculate the statistical p-value to achieve genome-wide association analysis, specifically including steps S51 to S53.

[0174] S51: Set the sample genotype vector of a certain site to be tested as G=(G 1T ,G2,…,G n ) T , the corresponding residual is R=(R1,R2,…,R n ) T .

[0175] S52: Under the assumption of Hardy Weinberg equilibrium, the genotypes of each member follow a binomial distribution, that is, G i ~B(2,μ), let is an estimate of μ, calculated based on the genetic effect β G The score statistic S under the assumption that is zero:

[0176]

[0177] S53: If The normal distribution approximation method is used to calculate the statistical p value; if The saddle point approximation method is used to calculate the statistical p-value, where is the observed value of S, is the estimated value of S variance. Step S53 is described in detail below.

[0178] (A1) Normal distribution approximation method

[0179] Under the assumption of Hardy-Weinberg equilibrium, the genotype variables of each member follow a binomial distribution, and the expectation and variance of S are:

[0180]

[0181]

[0182]

[0183] in, is the observed value of S, Φ(·) is the cumulative distribution function of the standard normal distribution, Var(G i ) represents the genotype variable G of the i-th individual i The variance of Cov(G i ,G j ) represents the genotype variable G of the i-th and j-th individualsi and G j The covariance of g Represents the genotype variable G i The standard deviation of μ represents the minor allele frequency of the tested locus, and ρ ij represents the kinship coefficient between members i and j, R i and R j Denotes the corresponding residual; let As an estimate of μ, As σ g Estimates, As an estimate of Var(S);

[0184] The expectation and variance of S are rewritten in matrix form as:

[0185] E(S)=2·μ·R T 1 n

[0186]

[0187] Among them, 1 n is an n-dimensional vector whose elements are all 1;

[0188] The statistical p value calculated using the normal distribution approximation method is:

[0189]

[0190] in, is the observed value of S, Φ(·) is the cumulative distribution function of the standard normal distribution, As an estimate of Var(S).

[0191] (A2) Saddle point approximation method

[0192] Unbalanced phenotypic and / or genotypic distributions can cause the normal distribution approximation of the score statistic to be inaccurate, significantly increasing the false positive rate. Using the saddlepoint approximation (SPA) method, we can obtain a more accurate estimate of the score statistic's distribution by leveraging information from higher-order moments, as follows.

[0193] Decompose the score statistic according to the sample's kinship coefficient and the corresponding residual:

[0194] S=S outliers +S non-outliers

[0195]

[0196] Among them, S outliersrepresents the score statistic of the family containing the sample data corresponding to the outlier, and there are K such families in total; S non-outliers The score statistic represents the family's sample data that does not contain outliers. The family information is not counted here. represents the genotype of the jth member of the i-th family containing the sample data corresponding to the outlier.

[0197] Using the saddle point approximation method requires estimating the moment generating function (MGF) of the score statistic. According to the above decomposition, the MGF of the score statistic is M S (t):

[0198]

[0199]

[0200] in, Indicates S non-putliers The moment generating function of Indicates S outliers The moment generating function of is approximated using the normal distribution method:

[0201]

[0202]

[0203]

[0204] Among them, E(S non-outliers ) indicates S non-outliers Expectations, Var(S non-outliers ) indicates S non-outliers The variance, R non-outliers represents the residual corresponding to the members of the family excluding the sample data corresponding to the outlier, Represents R non-outliers where μ represents the minor allele frequency.

[0205] Obtain the MGF of S (corresponding to m S (t)), the corresponding cumulative generating function (CGF) is:

[0206]

[0207] K S The first and second order derivative functions of (t) are:

[0208]

[0209]

[0210] Observed value of a given score statistic Compute the saddle point ζ to satisfy calculate and According to the Barndorff-Nielsen saddle point approximation formula, in Φ(·) represents the cumulative distribution function of the standard normal distribution, and then the two-sided statistical p value is obtained

[0211]

[0212] in yes Estimates, yes Estimates.

[0213] It should be noted that Stands for S outliers The calculation of MGF varies slightly depending on the number of family members, as follows:

[0214] If the family i containing the sample data corresponding to the outlier has only one member i1, and the genotype variable of member i1 follows the binomial distribution, we get in, represents the residual corresponding to member i1, and μ represents the minor allele frequency of the site to be tested.

[0215] If the family i containing the sample data corresponding to the outlier contains two members i1 and i2, the genotype variables of members i1 and i2 follow the binomial distribution, and the homologous sharing probability estimate of members i1 and i2 sharing 0, 1, and 2 genetic variations at the tested site is and the corresponding S i Moment generating function

[0216]

[0217]

[0218]

[0219] The sum of the products gives in Respectively represent the residuals corresponding to members i1, i2, S represents the case where members i1 and i2 share 0, 1, and 2 genetic variations at the same site, respectively. iThe moment generating function of , μ represents the minor allele frequency of the site to be tested;

[0220] If the family with the outlier corresponding to the sample data of the i-th family contains more than two members According to the estimated value of the joint probability of the genotype vectors of family members at the endpoints of each minor allele frequency interval, when the minor allele frequency μ of the site to be tested falls into a certain minor allele frequency interval, the joint probability distribution of the genotype vectors of the family members is It is approximated by the linear weighted average of the estimated values ​​of the joint probability at the two end points of the interval, S i Moment generating function

[0221]

[0222] By comparing the normal distribution approximation method and the saddle point approximation method described in (A1) and (A2), and One of the two methods is used to calculate the statistical p-value, realizing a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method.

[0223] According to steps S1 to S5, this analysis method fits a fitting constraint model that does not include genotypes when applied, that is, fits the constraint model and model parameters under the assumption that the genetic effect is zero and calculates the residuals; divides the residuals of the fitting constraint model into outliers and non-outliers, and divides the sample into K groups based on the kinship coefficient of the sample data. pre The invention provides a method for analyzing the genetic association of a whole genome in a fast and accurate manner. The invention provides a method for analyzing the genetic association of a whole genome in a fast and accurate manner. The method comprises the following steps: first, dividing the families with a number of members greater than c and containing sample data corresponding to the outliers into K families containing sample data corresponding to the outliers; second, calculating and storing the homology sharing probability between two members in the families with a number of members greater than 1 and containing sample data corresponding to the outliers; third, dividing all sites in the whole genome into a certain number of intervals according to the minor allele frequency (MAF); fourth, calculating and storing the joint probability distribution of the genotype vectors of each family at the endpoints of each interval using the Chow-Liu algorithm; fifth, estimating the MAF of each site to be tested and calculating the score statistic; fourth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; fifth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; fifth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; fifth, calculating the statistical p-value using a hybrid test strategy combining the normal distribution approximation method and the empirical saddle point approximation method to test the significance of the genetic effect; sixth ...

[0224] The following is a specific numerical simulation example:

[0225] This example provides a numerical simulation to demonstrate the effectiveness of the proposed method. This example applies the proposed algorithm to the analysis of kinship correlations in a control sample of longitudinal data, and numerical simulations are performed to evaluate the proposed method's Type I error rate and statistical power. The specific implementation steps are as follows:

[0226] (1) Generation of longitudinal data phenotypes

[0227] In the numerical simulation, three populations with a sample size of n = 50,000 are considered: small-family-based population A, consisting of 25,000 individuals and 6,250 families, each with 4 members; large-family-based population B, consisting of 25,000 individuals and 2,500 families, each with 10 members; and unrelated population C. The genetic pedigrees of the two families are shown in Figure 1. Figure 3 shown.

[0228] Longitudinal data phenotypes were generated using mixed-effects multi-location scale models:

[0229]

[0230]

[0231] Where 1≤i≤n. i Represents the longitudinal phenotypic observation of each individual, m i ~U(6,15). x ij and w ij Each contains three confounding factors: the first confounding factor does not change with time and is generated by Bernoulli distribution Bernoulli (0.5); the second confounding factor does not change with time and is generated by standard normal distribution; the third random variable changes with time and is generated by standard normal distribution. ij Contains a time-varying confounder generated by a standard normal distribution. i , W i and Z i All contain an intercept term with all values ​​1. The random effect term (γ i T ,ω i ) T Obey a mean vector of 0 q+1 , the covariance matrix is The multivariate normal distribution of random effect term b, Each obeys a mean vector of 0 n , the covariance matrix is ​​σ2 Sum of the multivariate normal distribution, where ∑ is an n×n dimensional genomic relationship matrix (GRM), and the GRM of the corresponding population is used.

[0232] Set β=(1, 0.5, 0.5, -0.3) T , τ=(0.25, 0.3, -0.15, 0.1) T , <00010​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​To evaluate the statistical power, other parameters remain unchanged and β G ≠0&τ G ≠ 0. Randomly select 5 common variant sites and 5 rare variant sites from the above 100,000 common variant sites and 100,000 rare variant sites, and let

[0239]

[0240] where g ik and For the kth normalized common variant and rare variant site, set the genetic effect to θ ik =-log10(MAF i )×0.08, and τ g =1.5,β g = 1. Each population structure produces 200 phenotypes, so 200 × 5 = 1000 tests of genetic effects are performed under each population structure and each variant type.

[0241] At the significance level of α = 5 × 10 -8 The statistical power of two methods is compared: a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method (SPA GRM ) and TrajGWAS methods. GRM ) In some cases, the first type error rate cannot be controlled, so the Norm GRM The TrajGWAS method cannot control the type I error rate in samples with kinship, so 25,000 independent samples from population A and 2 independent parents in each family, totaling 37,500 independent samples; 25,000 independent samples from population B and 4 independent parents in each family, totaling 35,000 independent samples; and all samples from population C were used for the TrajGWAS analysis. GRM All samples were analyzed acutely to compare the statistical power of the two methods.

[0242] (IV) Numerical simulation results of the type I error rate for longitudinal data phenotyping analysis:

[0243] The above is based on 10 8 The results of the numerical simulation of the first type error rate of the test are as follows Figure 4 As shown. Only the hybrid test strategy (SPA) combining the normal distribution approximation method and the saddle point approximation method GRM ) in any case, the first type error rate is controlled. G= 0, the TrajGWAS method can control the first type error rate in independent samples (i.e., population C), but cannot control the first type error rate in populations with kinship (i.e., populations A and B); G = 0, the TrajGWAS method can control the first type error rate when testing common variant sites, but it produces a large number of false positive sites when testing rare variant sites. GRM Only in the test β G = 0 and for common variant sites to control the first type error rate, while the other cases showed different degrees of overestimation. The results showed that under different kinship correlations and minor allele frequencies, SPA GRM The first type error rate can be controlled.

[0244] (V) Numerical simulation results of statistical power of longitudinal data phenotyping analysis:

[0245] The results of the numerical simulation of statistical power based on 1000 tests are as follows: Figure 5 As shown. At the significance level α=5×10 -8 Under these conditions, the two methods test β in an independent sample (i.e., population C). G = 0, the statistical power is almost the same when testing τ G =0 time SPA GRM The statistical power of the method is higher than that of the TrajGWAS method. In samples with kinship (i.e., population A and population B), since TrajGWAS can only analyze independent samples, it is not very effective in testing β G =0 and τ G =0 is less statistically powerful than SPA GRM .

[0246] The above numerical simulation experimental results show that the proposed hybrid test method combining normal distribution approximation and saddle point approximation (SPA GRM ) can well control the first type error rate; it can be applied to analyze samples with kinship to obtain higher statistical power.

[0247] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. An analytical method capable of controlling kinship correlation in a sample, characterized in that: The steps include: S1: Obtain sample data, wherein the sample data includes phenotypic data, genotypic data, and confounding factor data; S2: constructing a generalized linear model based on the sample data, ignoring the random effect terms in the generalized linear model, fitting the generalized linear model under the assumption that the genetic effect is zero to obtain a fitted constrained model, estimating the parameters of the fitted constrained model and calculating the residuals; S3: Divide the residuals of the fitted constraint model into outliers and non-outliers, and divide the sample data into K pre For each family, the number of family members is greater than c and the family contains sample data corresponding to the outlier is divided into K families containing sample data corresponding to the outlier. The probability of homology sharing between two members in the family with more than 1 member and containing sample data corresponding to the outlier is calculated and stored. S4: Divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency, and calculate and store the joint probability of the genotype vectors of each family member at the endpoints of each interval based on the homologous sharing probability and the Chow-Liu algorithm; S5: For each site to be tested, estimate the minor allele frequency and calculate the score statistic. Use a hybrid test strategy that combines the normal distribution approximation method with the empirical saddle point approximation method to calculate the statistical p-value to achieve genome-wide association analysis.

2. The method for analyzing kinship correlation in controllable samples according to claim 1, characterized in that: In step S2, it specifically includes: Assume that the number of sample data is n, let X represent an n×p dimensional confounding factor matrix, G represent the genotype vector, and Y represent a phenotypic data; Calculate the linear prediction term η=Xβ based on X, G, and Y X +Gβ G +b, treat the linear prediction term as a generalized linear model, where β X represents the p-dimensional coefficient vector of the confounding factor X, β G represents the genetic effect of genotype G, b represents the random effect term used to describe the relationship between individual members, b~N(0,τ∑), ∑ is the n×n dimensional genetic correlation matrix, τ is the additional genetic variance; Ignoring the random effect term b, the genetic effect β G Under the assumption that is zero, a generalized linear model is fitted to obtain a fitted constrained model, and the fitted constrained model is used to obtain the estimated values ​​of the parameters under the constrained conditions and to calculate the residuals of the fitted constrained model; Based on the estimated value, we can calculate the G The statistic S is zero.

3. The method for analyzing kinship correlation in controllable samples according to claim 2, characterized in that: In step S3, it specifically includes: Preselect a positive integer γ and let α 0.25 and α 0.75 They represent the 25% and 75% quantiles of the residual R obtained by fitting the constrained model, and the interquartile range is defined as IQR = α 0.75 -α 0.25 ; R i >α 0.75 +IRW*γ or R i <α 0.25 The residual corresponding to -IRQ·γ is regarded as an outlier, and the other residuals are regarded as non-outliers; The samples are divided into K groups based on the affinity coefficient of the sample data. pre A family, so that the kinship coefficient between family members is not 0, and the kinship coefficient between different family members is equal to 0; Preselect a positive integer c, split the families with more than c members and containing sample data corresponding to outliers, and obtain K families containing sample data corresponding to outliers; Based on any K families whose number of members is greater than 1 and contains the sample data corresponding to the outlier, set the i-th family to contain members 1,...,K i ,1≤i≤K, the corresponding residual is make represents the homologous sharing probability that the jth and kth members in the i-th family share 0, 1, and 2 genetic variants at the same locus, respectively; Under the assumption of Hardy-Weinberg equilibrium, the genotype variables of each member follow a binomial distribution; Using the information of s sites in the whole genome, we can obtain the kinship coefficient ρ for the jth and kth members. jk The probability of sharing 0 genetic variation at the same site with members j and k is λ jk Moment estimate of Moment-based estimation Calculate the probability of homology sharing Moment estimate of The moment estimate and moment estimates The calculation formula is as follows: in, represents the estimated value of the kinship coefficient between members j and k, represents the estimated probability of homologous sharing between members j and k sharing 0 genetic variants at the same site, G jl and G kl denote the genotype of the lth site of the jth and kth members, respectively. represents the weight of the genotype at the lth site, μ l represents the minor allele frequency of the lth site, represents μ l Estimates.

4. The method for analyzing kinship correlation in controllable samples according to claim 3, characterized in that: A positive integer c is pre-selected, and families with more than c members and containing sample data corresponding to outliers are segmented to obtain K families containing sample data corresponding to outliers, specifically including: Preselect a positive integer c; Assume that any family i with more than c members and containing outliers contains members The corresponding residual is Calculate the interaction parameter between any two members where ρ jk represents the kinship coefficient between members j and k, R j and R k represents the corresponding residual; The interaction parameter |ρ jk R j R k Arrange them in ascending order, and delete the related member pairs in sequence until the number of family members after the split is no more than c, and obtain K families that contain sample data corresponding to the outliers.

5. The method for analyzing kinship correlation in controllable samples according to claim 3, characterized in that: In step S4, it specifically includes: Divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency; Based on the fact that any family member number is greater than 2 and contains the sample data corresponding to the outlier, it is assumed that any family i contains members 1,...,K i ,K i >2, let G i Represents the genotype vector corresponding to the family members Let Pr(G i ) represents the joint probability of the genotype vector of the family member; Based on the genotype variable G of the jth and kth members j and G k Mutual information I(G j ,G k ), minimize the joint probability distribution Pr(G i ) and the approximate joint probability distribution The relative entropy of Minimizing the KL divergence is equivalent to calculating the j ,G k ),1≤j <k≤K i Construct (K i -1) Maximum spanning tree of edges; Using Prim's algorithm according to (K i (K i -1) / 2) mutual information I(G j ,G k ),1≤j <k≤K i , construct (K i -1) Maximum spanning tree of edges; Based on the Chow-Liu algorithm, according to (K i -1) The maximum spanning tree of edges obtains the optimal approximate joint probability distribution Approximate joint probability distribution As the genotype vector G of the family member i The joint probability of Repeat the above operation at each interval endpoint, calculate the joint probability of the genotype vectors of other family members, and store all the joint probabilities to obtain the joint probability of the genotype vectors of each family member at each interval endpoint corresponding to all sites in the whole genome.

6. The method for analyzing kinship correlation in controllable samples according to claim 5, characterized in that: described I(G j ,G k ) is calculated as follows: Among them, m l Indicates from 1 to K i Rearrangement, 1≤l≤K i , m h(l) represents the rearrangement from 0 to l-1, Indicates that at a given mth h(l) Member's genotype Under the condition m l Member's genotype The conditional probability of I(G j ,G k ) represents the genotype variable G of the jth and kth members j and G k The mutual information of H(G l ) is the genotype variable G of the lth member l The information entropy of H(G i ) represents the genotype vector of the family member The joint entropy of g j ,g k ∈{0,1,2} represents the genotype variable G of the jth and kth members respectively j and G k The possible values ​​of P(g j ),P(g k ) represent the genotype variables G of the jth and kth members respectively j and G k The observed value is g j ,g k The probability, P(g j ,g k ) represents the genotype variable G of the jth and kth members j and G k The observed value is g j ,g k The joint probability of .

7. The method for analyzing kinship correlation in controllable samples according to claim 3, characterized in that: In step S5, it specifically includes: Assume that the sample genotype vector of a certain site to be tested is G=(G1,G2,…,G n ) T , the corresponding residual is R=(R1,R2,…,R n ) T ; Under the assumption of Hardy-Weinberg equilibrium, the genotype variables of each member follow a binomial distribution, and the genetic effect β is calculated. G The score statistic S under the assumption that it is zero; (a) When using the normal distribution approximation method, the details are as follows: Calculate the expectation and variance of the score statistic S, and calculate the statistical p-value based on the expectation and variance: in, is the observed value of S, Φ(·) is the cumulative distribution function of the standard normal distribution, Var(G i ) represents the genotype variable G of the i-th individual i The variance of Cov(G i ,G j ) represents the genotype variable G of the i-th and j-th individuals i and G j The covariance of g Represents the genotype variable G i The standard deviation of μ represents the minor allele frequency of the site to be tested, and ρ ij represents the kinship coefficient between members i and j, R i and R j Denotes the corresponding residual; let As an estimate of μ, As σ g Estimates, As an estimate of Var(S); (b) Use the saddle point approximation method as follows: Decompose the score statistic S according to the kinship coefficient of the sample data and the corresponding residual to obtain the score statistic S of the family corresponding to the sample data of the outlier residual. outliers and the score statistic S of the family corresponding to the sample data without outlier residuals non-outliers ; According to S outliers and S non-outliers Calculate the moment generating function M of the score statistic S (t); The moment generating function M S (t) After accumulation, the moment function K is generated S (t); Based on the moment function K S (t) Calculate the statistical p value; in: S=S outliers +S non-outliers E(S non-outliers )2·μ·∑R non-outliers Among them, S i represents the score statistic of the family containing the sample data corresponding to the outlier, represents the genotype of the jth member of the i-th family containing the sample data corresponding to the outlier, Indicates S non-outliers The moment generating function of Indicates S outliers The moment generating function of represents the moment generating function of the score statistic of the family with sample data corresponding to the outlier, E(S non-outliers ) indicates S non-outliers Expectations, R non-outliers represents the residual corresponding to the members of the family excluding the sample data corresponding to the outlier, Represents R non-outliers The transpose of μ represents the minor allele frequency, Var(S non-outliers ) indicates S non-outliers The variance of represents the variance of the genotype variable at the site to be tested, ∑ non-outliers Represents the genetic correlation matrix between members of a family excluding sample data corresponding to outliers, where the jth and kth elements represent the kinship coefficient ρ of the jth and kth members jk , K′ S (t) and K″ S (t) represents the moment function K S The first and second order derivatives of (t), the saddle point ζ satisfies in, represents the observed value of S. According to the Barndorff-Nielsen saddle point approximation formula, in Φ(·) represents the cumulative distribution function of the standard normal distribution, and then the two-sided statistical p value is obtained in yes Estimates, yes estimates; (c) A hybrid testing strategy combining the normal distribution approximation method and the saddle point approximation method is used, as follows: Let r be a pre-selected positive number; if The normal distribution approximation method is used to calculate the statistical p value; if The saddle point approximation method is used to calculate the statistical p-value, where is the observed value of S, is the estimated value of the variance of S; After obtaining the statistical p-value, genome-wide association analysis was performed at a given significance level.

8. The method for analyzing kinship correlation in controllable samples according to claim 7, characterized in that: The calculations vary slightly depending on the number of family members, including: If the family i containing the sample data corresponding to the outlier has only one member i1, and the genotype variable of member i1 follows the binomial distribution, we get in, represents the residual corresponding to member i1, μ represents the minor allele frequency of the site to be tested; If the family i containing the sample data corresponding to the outlier contains two members i1 and i2, the genotype variables of members i1 and i2 follow the binomial distribution, and the homologous sharing probability estimate of members i1 and i2 sharing 0, 1, and 2 genetic variations at the tested site is and the corresponding S i Moment generating function of : The sum of the products gives in Respectively represent the residuals corresponding to members i1, i2, S represents the case where members i1 and i2 share 0, 1, and 2 genetic variations at the same site, respectively. i The moment generating function of , μ represents the minor allele frequency of the site to be tested; If the family with the outlier corresponding to the sample data of the i-th family contains more than two members According to the estimated value of the joint probability of the genotype variables of family members at the endpoints of each minor allele frequency interval, when the minor allele frequency μ of the site to be tested falls into a certain minor allele frequency interval, the joint probability distribution of the genotype variables of the family members is It is approximated by the linear weighted average of the estimated values ​​of the joint probability at the two end points of the interval, S i Moment generating function 9. An analysis system capable of controlling kinship correlation in a sample, characterized in that: It includes an acquisition module, a construction calculation module, a homology sharing probability calculation module, a joint probability calculation module and a statistical value calculation module; The acquisition module is used to acquire sample data, and the sample data includes phenotypic data, genotypic data and confounding factor data; The construction calculation module is used to construct a generalized linear model based on the sample data, ignore the random effect terms in the generalized linear model, fit the generalized linear model under the assumption that the genetic effect is zero to obtain a fitted constraint model, estimate the parameters of the fitted constraint model and calculate the residual; The homology sharing probability calculation module is used to divide the residuals of the fitted constraint model into outliers and non-outliers, and divide the sample data into K and N groups based on the kinship coefficient of the sample data. pre For each family, the number of family members greater than c and containing sample data corresponding to the outlier is divided into K families containing sample data corresponding to the outlier. The probability of homology sharing between two members in the family with more than 1 and containing sample data corresponding to the outlier is calculated and stored. The joint probability calculation module is used to divide all sites in the whole genome into a certain number of intervals according to the minor allele frequency, and calculate and store the joint probability of the genotype vectors of each family member at the endpoints of each interval based on the homologous sharing probability and the Chow-Liu algorithm; The statistical value calculation module is used to estimate the minor allele frequency and calculate the score statistic for each site to be tested, and calculate the statistical p-value using a hybrid test strategy that combines the normal distribution approximation method with the empirical saddle point approximation method to achieve genome-wide association analysis.

10. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a plurality of programs, which are used to be called by a processor and execute the method for analyzing kinship correlation in controllable samples as claimed in any one of claims 1 to 8.

Images