Gene level whole genome association analysis algorithm based on EMS population

Abstract

The invention discloses a gene level whole genome association analysis algorithm based on an EMS population, and relates to the technical field of bioinformatics. The correlation analysis is carried out by taking the gene as a basic unit, the MAF index of the SNP of the EMS population is far lower than that of a normal population, but the MAF value of the gene level is higher than that of the SNP level. According to the method, weighting is carried out according to the mutation effect of each mutation site, and the total weighted value of all mutations in each gene of a single sample is counted as the basis of correlation analysis, so that the false positive rate of the result is greatly reduced, and the statistical effectiveness of analysis is improved. Meanwhile, a plurality of statistical methods are used for comprehensive judgment so as to find the most reliable candidate gene. Compared with the situation that GWAS can only locate one fuzzy interval possibly containing multiple genes, the method can directly and accurately locate a single gene, the experimental efficiency of candidate gene function verification is improved, and the problems that in the prior art, large-scale analysis of an EMS population cannot be borne, and the analysis efficiency is low are effectively solved.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a gene-level genome-wide association analysis algorithm based on an EMS population. Background Art

[0002] EMS mutagenesis (Ethyl Methanesulfonate) is a commonly used chemical mutagenesis method used to cause base mutations in biological DNA through chemical substances. The EMS mutant library is a large collection of mutants generated by using EMS to mutate organisms (such as plants, animals, bacteria, etc.). The characteristics of the EMS mutant library are high throughput and wide distribution of mutations, but the mutation positions caused by EMS mutagenesis are usually random. This results in the Minor Allele Frequency (MAF) of each base site not meeting the requirements for genome-wide association studies (GWAS), which brings difficulties to GWAS analysis. GWAS is to find SNPs across the entire genome, screen out SNPs related to traits, and study the correlation between genetic mutations and phenotypes. It can be used to study important traits such as crop yield and disease resistance, providing an important basis for agricultural breeding.

[0003] In human disease research, GWAS analysis mainly uses large sample sizes to find differences in common variants (MAF>1%) between case and control groups to explain complex diseases, ignoring the role and information of rare variants (MAF<1%), which are believed to explain the genetic relationship and mechanism of some diseases. Therefore, a region / gene-based association analysis method, the Burden analysis, was developed. Simply put, the Burden test compares two groups of samples with different phenotypes to see if there is a significant difference in the "total number" of rare variants carried in the same region / gene.

[0004] However, the load test cannot be applied to the EMS mutant library. First, the number of rare variants in humans is much lower than the number of EMS sites, with a difference of millions. The thousands of individuals in the EMS mutant library are also much larger than the hundreds of cases in the human cohort, and the load test cannot afford large-scale analysis. Secondly, the phenotype of the disease is a qualitative trait, while the phenotype of crops is a quantitative trait. Using the load test to analyze quantitative traits will greatly reduce the effectiveness of the analysis. Thirdly, the load test counts the number of mutations in the statistical area, and there is no reasonable weighting of mutation sites with different functions, resulting in reduced analysis effectiveness. Finally, the load test also requires a large number of control samples, which increases the threshold for analysis. Summary of the invention

[0005] In view of the above-mentioned deficiencies in the prior art, the present invention provides a gene-level whole-genome association analysis algorithm based on the EMS population. The present invention can directly and accurately locate a single gene, improve the efficiency of the candidate gene function verification experiment, and effectively solve the problems that the prior art cannot afford large-scale analysis, has low analysis efficiency and high analysis threshold.

[0006] To achieve the above object, the technical solution adopted by the present invention to solve the technical problem is: to provide a gene-level genome-wide association analysis algorithm based on the EMS population, comprising the following steps: S1. Obtain the sample to be analyzed for exome sequencing; S2, obtain phenotypic data of EMS samples; S3, compare the exome sequencing data of the EMS population with the genome of its species; S4, using a SNP detection tool to perform SNP detection on the obtained comparison results to obtain SNP information in the gene coding region; S5. Use SnpEff software to annotate SNPs and predict mutation effects to obtain the degree of influence of each SNP variation on gene function; S6. Count the total mutation impact on each gene region of each sample to obtain the impact degree information at the gene level; S7. Association analysis between genotype and phenotype was performed using linear models and chi-square test.

[0007] Furthermore, in step S3, the exon sequencing data in step S1 are compared with the whole genome of the species using comparison software such as BWA and Bowtie2.

[0008] Furthermore, in step S4, SNP detection is performed on the obtained alignment results using SNP detection tools such as GATK, Samtools, BCFtools, etc.

[0009] Further, in step S4, the SNP information is saved in a VCF file format.

[0010] Further, in step S4, the specific method is: extract the genotype information of the SNP sites of all samples from the VCF file, obtain the typing matrix of the SNP sites, and filter out abnormal sites with too high mutation frequency according to the mutation frequency threshold determined by the permutation test.

[0011] Furthermore, the genotype information and sample information of each SNP in the VCF file are shuffled to obtain randomly arranged sample genotype information, and the mutation frequency of each site is counted and arranged from low to high. The 5th lowest frequency value is selected as the result of a test. This is performed n times in total, and the n results are arranged from low to high, and the n*0.01 lowest result is selected as the mutation frequency threshold.

[0012] Further, in step S5, the specific method is: use the annotation tool SnpEff to annotate the specific gene where each site is located, set the distance threshold K, K is set to a positive number less than 3000, obtain the SNP sites on the gene and within K bp upstream and downstream of the gene as the SNP sites in the gene region, and predict the impact effect of each SNP, and the levels of the impact effect are divided into HIGH, MODERATE, LOW and MODIFIER.

[0013] Further, in step S7, a genotype-phenotype association analysis is performed, which specifically includes the following steps: S71, filtering false positive SNPs; S72, filter phenotypic outliers; S73. Establish a positive statistical model of gene-level mutation effect and phenotype, conduct association analysis between genotype and phenotype, and obtain the association P value between phenotype and gene; S74. Establish a reverse statistical model of the mutant phenotype and the mutation effect at the gene level, conduct association analysis between the phenotype and the genotype, and obtain the P value of the association between the gene and the phenotype; S75, merging the forward and reverse results to obtain the final association score between the phenotype and the gene; S76. Based on the principle of statistical testing, use permutation testing to obtain a reliable significance threshold; S77. Screening of significantly associated genes.

[0014] Further, in step S73, the association analysis between genotype and phenotype is performed, which specifically includes the following steps: S731. Each site is scored according to the level of the site mutation effect, where MODIFIER=1.000001, LOW=2.00001, MODERATE=3.0001, HIGH=4.001; S732. According to the genotype information and gene position annotation of the sample, the total mutation effect score of each gene of each sample is calculated. For example, if there are 3 LOW level mutations and 2 HIGH level mutations on gene A of a sample, the total mutation effect score of gene A is 14.00203. S733. Establish a linear regression statistical model of genotype and phenotype at the gene level, perform association analysis between genotype and phenotype, obtain the association P value between phenotype and each gene, and select the top 0.5% as significantly associated genes; S734. Establish a test model for gene mutation and phenotype, conduct association analysis between gene mutation and phenotypic abnormalities, obtain the P value of association between phenotype and each gene, and select the top 0.5% as significantly associated genes; S735, verifying the significant associated genes in step S33 and step S34, and merging and taking the intersection as candidate associated genes; S736. Determine the significance threshold through permutation test, filter the candidate associated genes whose significance is less than the threshold, and obtain the final result.

[0015] Further, in step S733, for a single gene A, the total mutation effect score of gene A in each sample is calculated as the linear model x, the phenotypic value is taken as y, and a linear regression fitting of y=ax+b is performed on the entire EMS population to determine the significance P value of the linear model fitting coefficient as the P value of the association between gene A and phenotype. Statistics are performed on all genes to obtain the final result.

[0016] Further, in step S734, for gene A, samples with significant phenotypic differences are determined by box plots, the number is counted, and the number of samples with mutations and the number of samples without mutations in gene A are counted. The significance relationship between the number of samples with abnormal phenotypes and the number of samples with mutations in gene A is obtained by chi-square test, which is used as the P value of the association between gene A and phenotype. Statistics are performed on all genes to obtain the final result.

[0017] Furthermore, in step S736, after the phenotypic information of the population and the sample information are disrupted, the association analysis of S734 and S735 is performed again, and the P value of the fifth gene with the highest significance is recorded. This is repeated n times, and the n results are arranged from low to high, and the result with the lowest n*0.01 is selected as the significance threshold.

[0018] Further, in step S77, the significantly associated genes are screened, specifically comprising the following steps: S771, setting a distance threshold K, where K is set to a positive number less than 3000, and obtaining SNP sites on the gene and within K bp upstream and downstream of the gene; S772, using simulation data to predict the mutation frequency threshold of SNPs in the population, and filtering SNP sites that exceed the threshold; S773, performing comprehensive calculation based on the gene association P values ​​calculated in step S73 and step S74 to obtain a final gene association score; S774. Determine a threshold p0 of the association score, and screen out genes with an association score greater than p0.

[0019] The present invention has the following beneficial effects: 1. The genotype of the EMS mutation population has the characteristics of sparseness and randomness. The sparseness makes the MAF of a single site too low, so that it is impossible to use the existing association analysis technology to perform genotype-phenotype association analysis. The basic unit of the genotype of the association analysis is increased from SNP to the gene level, which greatly improves the effectiveness of the detection. At the same time, the present invention is not a simple statistical gene mutation, but a score for the overall effect of the gene according to the mutation effect of the gene, so as to construct a linear regression model and improve the accuracy of the detection.

[0020] 2. The present invention also uses mutation enrichment test to reversely verify the results, reducing the error rate of the results. According to the complexity of data distribution and model, the traditional critical value is not used as the threshold, but the permutation test result is used as the threshold to prevent false positive results. When performing association analysis at the gene level, the results obtained are directly specific genes, rather than the genomic regions of traditional association analysis, which greatly reduces manpower, material resources and time costs.

[0021] 3. Association analysis is carried out with genes as basic units, and the MAF index is greatly improved compared to the base level. Weighted according to the effect of each mutation site, the total weighted value of all mutations in the single gene of a single sample is counted, as the basis of association analysis, the false positive rate of the result is greatly reduced, and the statistical effectiveness of the analysis is improved. A variety of statistical methods are used for comprehensive evaluation at the same time, to find the gene most relevant to the phenotype. At the same time, relative to GWAS, only a vague interval can be located, and the result may include multiple genes. The present invention can directly and accurately locate a single gene, which improves the efficiency of the candidate gene function verification experiment. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 Schematic diagram of determining filtering thresholds for permutation tests; Figure 2 Schematic diagram of phenotypic outliers and phenotypic mutations; Figure 3 Schematic diagram of non-significant genes obtained by linear model association; Figure 4 Schematic diagram of significant genes obtained by linear model association. DETAILED DESCRIPTION

[0023] The principles and features of the present invention are described below. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0024] Example 1 The EMS mutation population data of wheat variety KN9204 are used to illustrate the genotype and phenotype association analysis method based on EMS population described in the present invention.

[0025] S1. 2000 seeds of wheat KN9204 selected in our laboratory in 2020 were induced with EMS and planted in the field. 1800 strains of wheat KN9204 induced population that survived in the field were sampled for single-plant whole exon sequencing. Whole exon capture sequencing was performed using a coding sequence panel designed for the KN9204 genome, including 1232636 probes covering 127 Mb of the KN9204 genome interval. The sequencing platform was BGI T7, and the sequencing depth was 30X. S2, harvesting the surviving strains described in step S1 and performing phenotypic identification, measuring the thousand-grain weight of each strain, and obtaining phenotypic data of all samples, the results are shown in Table 1; Table 1 1810 sample thousand-grain weight data

[0026] S3. Use BWA software to align the whole exome sequencing data of 1800 samples (default parameters), and the reference genome is the published KN9204 genome; S4. The GATK tool (https: / / gatk.broadinstitute.org / ) was used to perform SNP detection on the parental alignment results. A total of 5.9M sites were detected. The BCFtools tool (https: / / samtools.github.io / bcftools / ) was used for screening (the screening index was %QUAL<50 || INFO / DP<5). The permutation test was used to determine the threshold of mutation frequency as 0.007, and false positive sites with too high mutation frequency were filtered out, such as Figure 1 As shown; The genotype information of the SNP sites of the samples was extracted from the VCF file (for example, the sample numbered M1800 was extracted to have a heterozygous mutant genotype at Chr1A_3920, and the sample numbered M0002 was extracted to have a homozygous mutant genotype at Chr7D_615206503), and the genotype matrix of the SNP sites in the samples was obtained. The results are shown in Table 2; Table 2 Genotype data of 1800 samples

[0027] S5. Use SnpEff software to annotate SNPs and predict mutation effects, obtain the degree of influence of each SNP variation on gene function, use SnpEff tool to analyze the location of each site on the genome, obtain the gene and mutation effect level of each site, and assign points according to the effect level. The results are shown in Table 3. Table 3 Mutation effect of 5.9M SNP

[0028] S6. Count the total mutation impact on each gene region of each sample to obtain the impact degree information at the gene level; S7. Association analysis between genotype and phenotype was performed using linear models and chi-square test.

[0029] The above-mentioned method for association analysis between genotype and phenotype comprises the following steps: S71. Filter out false positive sites with too high mutation frequency in the VCF file. After randomly matching the genotype information of each SNP in the genotype data with the sample information, count the mutation frequency of each site, sort them from low to high, and select the fifth result as the result of the current test. Perform a total of 1000 times, sort the 1000 results from low to high, and select the 10th result as the mutation frequency threshold. SNP sites that exceed the mutation frequency threshold are filtered out. S72. Filtering of outlier phenotypes in mutant populations: Use R software to draw box plots, such as Figure 2 As shown, 2 extreme outliers were identified, and the phenotypic and genotypic data of the remaining 1798 EMS population samples were retained for subsequent analysis; S73. Calculate the overall score of the mutation effect of each sample in each gene based on the genotype annotation information to obtain the sample genotype mutation effect score matrix. The results are shown in Table 4. Table 4 Sample genotype mutation effect scores

[0030] Among them, samples without mutation sites on the gene had the mutation effect of the gene recorded as 0 and were included in the analysis.

[0031] S74. Based on the generalized linear model, the phenotypic data and gene mutation effect score data of each sample were combined to perform an association analysis between the genotype mutation effect score and its thousand-grain weight in the R environment, such as Figure 3 and Figure 4 As shown, with the mutation effect score as x and the thousand-grain weight data as y, a generalized linear model y=ax+b was constructed, and the model was regressed to obtain the fitted linear model coefficients. The significance level of the coefficients was obtained through the t-test, and the associated P value of the gene was obtained; S75, use Figure 2The box plot drawn by R software shown in the figure determines 65 samples with significant differences outside the interval [Q1-3IQR, Q3+3IQR], where Q1 and Q3 represent the lower quartile and upper quartile of the sample thousand-grain weight data, respectively, and IQR represents the interquartile range. All samples are grouped according to the mutation effect score of the gene, one group with a mutation effect score of 0 and the other group with a mutation effect score of not 0. If the gene is associated with the phenotype, the phenotype of the sample in the group with a mutation effect score of not 0 should also be significantly different from the whole. The chi-square test is used to test whether the phenotype of the mutant samples of each gene is significantly different, and the associated P value of the gene is obtained; S76. Perform a permutation test to calculate the significance threshold. After disrupting the correspondence between the sample phenotype and the sample number, perform steps S73 and S74 to calculate, respectively count the value with the smallest P value, repeat 1000 times, and obtain the 1000 minimum P value results of the two methods respectively. Sort them from small to large according to the P value, select the 50th ranked P value as the significance threshold, calculate the -log value of the average of the two p values, and finally determine the threshold to be 3.2.

[0032] Take TraesKN3D01HG16320 as an example to illustrate the results of the above steps: The situation of TraesKN3D01HG16320 is shown in Table 5.

[0033] Table 5 The situation of TraesKN3D01HG16320

[0034] The significance threshold in step S74 is 0.000089 and in step S74 is 0.000103.

[0035] The P value of association between TraesKN3D01HG16320 and the phenotype calculated using the method in step S73 was 0.000032, ranking 103rd among all genes (a total of 108,326 genes).

[0036] The P value of association between TraesKN3D01HG16320 and phenotype calculated using the method in step S74 is 0.000059, ranking 166th among all genes (a total of 108,326 genes). The average -log value of the P values ​​in the two steps is 4.34.

[0037] Therefore, this gene was identified as a candidate gene related to thousand-grain weight.

[0038] The list of highly supported group-associated genes based on the novel association analysis method is shown in Table 6.

[0039] Table 6 List of highly supported group-associated genes based on the new association analysis method

[0040] As shown in Table 6, thousand-grain weight is a complex trait affected by multiple factors. The genes in Table 6 achieve the regulation of wheat growth and development from multiple angles, which is consistent with the developmental characteristics of the population. The method provided by the present invention is used in various crop populations induced by EMS, such as rice, wheat, corn, etc. The target trait can be a quantitative trait or a multivariate quality trait.

[0041] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A gene-level genome-wide association analysis algorithm based on the EMS population, characterized in that: The following steps are involved: S1. Obtain the sample to be analyzed for exome sequencing; S2, obtaining phenotypic data of EMS samples and normal control samples; S3, compare the exome sequencing data of the EMS population with the genome of its species; S4, using a SNP detection tool to perform SNP detection on the obtained comparison results to obtain SNP information in the gene coding region; S5. Use SnpEff software to annotate SNPs and predict mutation effects to obtain the degree of influence of each SNP variation on gene function; S6. Count the total mutation impact on each gene region of each sample to obtain the impact degree information at the gene level; S7. Association analysis between genotype and phenotype was performed using linear models and chi-square test.

2. The gene-level genome-wide association analysis algorithm based on the EMS population according to claim 1, characterized in that: In step S4, the SNP information is saved in a VCF file format. The specific method is: extract the genotype information of the SNP sites of all samples from the VCF file to obtain the typing matrix of the SNP sites, and filter out abnormal sites with too high mutation frequency according to the mutation frequency threshold determined by the arrangement test.

3. The gene-level genome-wide association analysis algorithm based on the EMS population according to claim 1, characterized in that: In step S5, the specific method is: use the annotation tool SnpEff to annotate the specific gene where each site is located, set the distance threshold K, K is set to a positive number less than 3000, obtain the SNP sites on the gene and within K bp upstream and downstream of the gene as the SNP sites in the gene region, and predict the impact effect of each SNP. The levels of the impact effect are divided into HIGH, MODERATE, LOW and MODIFIER.

4. The EMS population-based gene-level genome-wide association analysis algorithm according to claim 1, characterized in that: In step S7, the association analysis between genotype and phenotype is performed, which specifically includes the following steps: S71, filtering false positive SNPs; S72, filter phenotypic outliers; S73. Establish a positive statistical model of gene-level mutation effect and phenotype, conduct association analysis between genotype and phenotype, and obtain the association P value between phenotype and gene; S74. Establish a reverse statistical model of the mutant phenotype and the mutation effect at the gene level, conduct association analysis between the phenotype and the genotype, and obtain the association P value between the gene and the phenotype; S75, merging the forward and reverse results to obtain the final association score between the phenotype and the gene; S76. Based on the principle of statistical testing, use permutation testing to obtain a reliable significance threshold; S77. Screening of significantly associated genes.

5. The gene-level genome-wide association analysis algorithm based on the EMS population according to claim 4, characterized in that: In step S73, the association analysis between genotype and phenotype is performed, which specifically includes the following steps: S731. Each site is scored according to the level of the site mutation effect, where MODIFIER=1.000001, LOW=2.00001, MODERATE=3.0001, HIGH=4.001; S732. According to the genotype information and gene position annotation of the sample, the total mutation effect score of each gene of each sample is calculated. For example, if there are 3 LOW level mutations and 2 HIGH level mutations on gene A of a sample, the total mutation effect score of gene A is 14.00203. S733. Establish a linear regression statistical model of genotype and phenotype at the gene level, perform association analysis between genotype and phenotype, obtain the association P value between phenotype and each gene, and select the top 0.5% as significantly associated genes; S734. Establish a test model for gene mutation and phenotype, conduct association analysis between gene mutation and phenotypic abnormalities, obtain the P value of association between phenotype and each gene, and select the top 0.5% as significantly associated genes; S735, verifying the significant associated genes in step S3-3 and step S3-4, and merging and taking the intersection as candidate associated genes; S736. Determine the significance threshold through permutation test, filter the candidate associated genes whose significance is less than the threshold, and obtain the final result.

6. The EMS population-based gene-level genome-wide association analysis algorithm according to claim 4, characterized in that: In step S77, the significantly associated genes are screened, specifically comprising the following steps: S771, setting a distance threshold K, where K is set to a positive number less than 5000, and obtaining SNP sites on the gene and within K bp upstream and downstream of the gene; S772, using simulation data to predict the mutation frequency threshold of SNPs in the population, and filtering SNP sites that exceed the threshold; S773, performing comprehensive calculation according to the gene association P values ​​calculated in step S73 and step S74 in claim 4 to obtain a final gene association score; S774. Determine a threshold p0 of the association score, and screen out genes with an association score greater than p0.

Images