Genome-wide Association Analysis Algorithm at the Gene Level Based on the EMS Population

Abstract

The present invention discloses a genome-wide association analysis algorithm at the gene level based on an EMS population, which relates to the technical field of bioinformatics. The present invention conducts association analysis with genes as the basic unit. The MAF index of SNPs in the EMS population is much lower than that of the normal population, but the MAF value at the gene level is higher than that at the SNP level. Weights are assigned according to the mutation effects of each mutation site, and the total weighted value of all mutations in each gene of a single sample is statistically calculated as the basis for association analysis, greatly reducing the false positive rate of the results and increasing the statistical power of the analysis. At the same time, multiple statistical methods are used for comprehensive evaluation to find the most reliable candidate genes. Compared with GWAS that can only locate a fuzzy interval that may contain multiple genes, the present invention can directly and accurately locate to a single gene, improving the experimental efficiency of functional verification of candidate genes and effectively solving the problems that the prior art cannot afford large-scale analysis of the EMS population and has low analysis power.

Description

Technical Field

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

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

[0003] In GWAS analysis for human disease research, it mainly uses a large sample size 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%). Rare variants are considered to be able to explain part of the genetic relationship and pathogenesis of diseases. Therefore, a burden analysis, a region / gene-based association analysis method, has been developed. Briefly speaking, the burden analysis is to compare whether there is a significant difference in the "total number" of rare variants carried in the same region / gene between two samples with different phenotypic differences.

[0004] However, the burden analysis 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 of levels. The number of individuals in the EMS mutant library is also much larger than the few hundred cases in the human cohort, and the burden analysis cannot handle large-scale analysis. Second, the phenotypes of diseases are qualitative traits, while the phenotypes of crops are quantitative traits. Using the burden analysis to analyze quantitative traits will greatly reduce the analysis power. Third, the burden analysis only counts the number of mutations in the statistical region and does not reasonably weight mutation sites with different functions, resulting in a decrease in analysis power. Finally, the burden analysis also requires a large number of control samples, increasing the analysis threshold. Summary of the Invention

[0005] In view of the above deficiencies in the prior art, the present invention provides a genome-wide association analysis algorithm at the gene level based on an EMS population. The present invention can directly and accurately locate a single gene, improving the efficiency of candidate gene function verification experiments and effectively solving the problems of the prior art that large-scale analysis cannot be afforded, the analysis power is low, and the analysis threshold is high.

[0006] To achieve the above object, the technical solution adopted by the present invention to solve its technical problems is: to provide a genome-wide association analysis algorithm at the gene level based on an EMS population, including the following steps:

[0007] S1. Obtain the exome sequencing of the sample to be analyzed;

[0008] S2. Obtain the phenotypic data of the EMS sample;

[0009] S3. Align the exome sequencing data of the EMS population with the genome of its species;

[0010] S4. Use SNP detection tools to perform SNP detection on the obtained alignment results to obtain SNP information in the gene coding region;

[0011] S5. Use SnpEff software to annotate SNPs and predict mutation effects to obtain the degree of influence of each SNP variation on gene function;

[0012] S6. Statistically analyze the total mutation effects received by each gene region of each sample to obtain information on the degree of influence at the gene level;

[0013] S7. Use linear models and chi-square tests for genotype and phenotype association analysis.

[0014] Further, in step S3, alignment software such as BWA and Bowtie2 is used to align the exome sequencing data in step S1 with the whole genome of its species.

[0015] Further, in step S4, SNP detection tools such as GATK, Samtools, and BCFtools are used to perform SNP detection on the obtained alignment results.

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

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

[0018] Further, shuffle the genotype information of each SNP in the VCF file with the sample information to obtain randomly arranged sample genotype information, count the mutation frequency at each locus, arrange them from low to high, select the fifth lowest frequency value as the result of one test, perform this n times, arrange the n results from low to high, and select the result at the n*0.01 lowest position as the mutation frequency threshold.

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

[0020] Further, in step S7, perform the association analysis between genotype and phenotype, which specifically includes the following steps:

[0021] S71. Filter false positive SNPs;

[0022] S72. Filter phenotypic outliers;

[0023] S73. Establish a positive statistical model between the gene-level mutation effect and the phenotype, perform the association analysis between genotype and phenotype, and obtain the association P value between the phenotype and the gene;

[0024] S74. Establish a reverse statistical model between the mutant phenotype and the gene-level mutation effect, perform the association analysis between the phenotype and the genotype, and obtain the association P value between the gene and the phenotype;

[0025] S75. Combine the forward and reverse results to obtain the final association score between the phenotype and the gene;

[0026] S76. According to the principle of statistical test, use permutation testing to obtain a reliable significance threshold;

[0027] S77. Screen significantly associated genes.

[0028] Further, in step S73, the association analysis between genotype and phenotype specifically includes the following steps:

[0029] S731. Assign scores to each locus according to the impact effect level of the locus mutation. Among them, MODIFIER = 1.000001, LOW = 2.00001, MODERATE = 3.0001, HIGH = 4.001;

[0030] S732. According to the genotype information and gene position annotation of the samples, calculate the total mutation impact effect score of each gene in each sample. For example, if there are 3 LOW-level mutations and 2 HIGH-level mutations on gene A of a sample, then the total mutation effect score of gene A is 14.00203;

[0031] S733. Establish a linear regression statistical model of genotype and phenotype at the gene level, conduct an association analysis of genotype and phenotype, obtain the association P-value between phenotype and each gene, and select the top 0.5% as significantly associated genes;

[0032] S734. Establish a test model for the gene mutation situation and the phenotype situation, conduct an association analysis of gene mutation and phenotype abnormality, obtain the association P-value between phenotype and each gene, and select the top 0.5% as significantly associated genes;

[0033] S735. Mutually verify the significantly associated genes in steps S33 and S34, and take the intersection after merging as candidate associated genes;

[0034] S736. Determine the significance threshold through permutation test, filter the candidate associated genes with significance less than the threshold, and obtain the final result.

[0035] Furthermore, in step S733, for a single gene A, calculate the total mutation effect score of gene A in each sample as the linear model x, and the phenotype value as y. Conduct a linear regression fitting of y = ax + b for 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 the phenotype. Conduct statistics on all genes to obtain the final result.

[0036] Furthermore, in step S734, for gene A, determine the samples with significant phenotypic differences through a box plot, count the number, and at the same time count the number of samples with mutations on gene A and the number of samples without mutations. Obtain the significant relationship between the number of samples with abnormal phenotypes and the number of samples with mutations on gene A through a chi-square test as the P-value of the association between gene A and the phenotype. Conduct statistics on all genes to obtain the final result.

[0037] Furthermore, in step S736, after shuffling the phenotypic information and sample information of the population, conduct the association analysis of S734 and S735, and record the P-value of the fifth gene with the highest significance. Repeat this n times, arrange the n results from low to high, and select the result of the n * 0.01 lowest as the significance threshold.

[0038] Furthermore, in step S77, screen the significantly associated genes, which specifically include the following steps:

[0039] S771. Set a distance threshold K, where K is set as a positive number less than 3000, and obtain the SNP sites on the gene and within K bp upstream and downstream of the gene;

[0040] S772. Use simulated data to predict the mutation frequency threshold of SNPs in the population, and filter out the SNP sites exceeding the threshold;

[0041] S773. Based on the gene association P values calculated in steps S73 and S74, perform comprehensive calculations to obtain the final gene association score;

[0042] S774. Determine the threshold p0 of the association score, and screen out the genes with an association score greater than p0.

[0043] The present invention has the following beneficial effects:

[0044] 1. The genotypes of the EMS mutant population have the characteristics of sparsity and randomness. Sparsity makes the MAF of a single site too low to perform genotype-phenotype association analysis using existing association analysis techniques. By raising the basic unit of the genotype for association analysis from SNP to the gene level, the detection effectiveness is greatly improved. At the same time, the present invention does not simply count whether there is a mutation in the gene, but scores the overall effect of the gene according to the mutation effect of the gene, thereby constructing a linear regression model and improving the detection accuracy.

[0045] 2. The present invention also uses mutation enrichment tests to reverse-verify the results, reducing the error rate of the results. According to the data distribution and the complexity of the model, instead of using the traditional critical value as the threshold, the permutation test results are used as the threshold to prevent false positives in the results. By performing association analysis at the gene level, the obtained results are directly specific genes, rather than genomic regions in traditional association analysis, greatly reducing the human, material and time costs.

[0046] 3. When performing association analysis with the gene as the basic unit, the MAF index is greatly improved compared to the base level. By weighting according to the effect of each mutation site and statistically calculating the total weighted value of all mutations in a single gene of a single sample as the basis for association analysis, the false positive rate of the results is greatly reduced, and the statistical power of the analysis is improved. At the same time, multiple statistical methods are used for comprehensive evaluation to find the gene most relevant to the phenotype. At the same time, compared with GWAS which can only locate to a vague interval and the results may contain multiple genes, the present invention can directly and accurately locate to a single gene, improving the efficiency of candidate gene function verification experiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 Schematic diagram for determining the filtering threshold by permutation test;

[0048] Figure 2Schematic diagram of phenotypic outliers and phenotypic mutations;

[0049] Figure 3 Schematic diagram of non-significant genes obtained by linear model association;

[0050] Figure 4 Schematic diagram of significant genes obtained by linear model association. Detailed implementation manners

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

[0052] Example 1

[0053] Taking the data of the EMS mutant population of wheat variety KN9204 as an example to illustrate the genotype and phenotype association analysis method based on the EMS population described in the present invention.

[0054] S1. 2000 seeds of individual plants of wheat KN9204 selected by this laboratory in 2020 were subjected to EMS mutagenesis and planted in the field; 1800 surviving wheat KN9204 mutant population lines in the field were sampled for whole-exome sequencing of individual lines. Whole-exome capture sequencing was performed using a coding sequence panel specifically designed for the KN9204 genome, including 1,232,636 probes covering a 127 Mb KN9204 genome interval. The sequencing platform was BGI T7, and the sequencing depth was 30X;

[0055] S2. After harvesting the surviving lines described in step S1, phenotypic identification was carried out, and the thousand-grain weight of each line was measured to obtain the phenotypic data of all samples. The results are shown in Table 1;

[0056] Table 1 Thousand-grain weight data of 1810 samples

[0057]

[0058] S3. The BWA software was used to align the whole-exome sequencing data of 1800 samples (default parameters), and the reference genome was the published KN9204 genome;

[0059] S4. The GATK tool (https: / / gatk.broadinstitute.org / ) was used to detect SNPs in the obtained parental alignment results. A total of 5.9M loci were detected. The BCFtools tool (https: / / samtools.github.io / bcftools / ) was used for screening (the screening index was %QUAL<50 || INFO / DP<5); The threshold of the mutation frequency was determined by permutation test to be 0.007, and false positive sites with too high mutation frequency were filtered, such asFigure 1 as shown

[0060] Extract the genotype information of SNP sites of samples from the VCF file (for example, it is extracted that the sample numbered M1800 has a heterozygous mutant genotype at Chr1A_3920, and the sample numbered M0002 has a homozygous mutant genotype at Chr7D_615206503), and obtain the genotype matrix of SNP sites in the samples. The results are shown in Table 2;

[0061] Table 2 Genotype data of 1800 samples

[0062]

[0063] S5. Use the SnpEff software to annotate SNPs and predict mutation effects, obtain the degree of influence of the variation of each SNP on gene function, use the SnpEff tool to analyze the positions of each site on the genome, obtain the gene where each site is located and its own mutation effect level, and assign scores according to the effect level. The results are shown in Table 3;

[0064] Table 3 Mutation effects of 5.9M SNPs

[0065]

[0066] S6. Statistically analyze the total mutation effects received by each gene region of each sample to obtain the information on the degree of influence at the gene level;

[0067] S7. Use linear models and chi-square tests for genotype-phenotype association analysis.

[0068] The method for genotype-phenotype association analysis described above includes the following steps:

[0069] S71. Filter out false positive sites with too high mutation frequencies in the VCF file. After randomly matching the genotype information of each SNP in the genotype data with the sample information, calculate the mutation frequency of each site, arrange them from low to high, select the fifth result as the result of the current test, and perform this 1000 times. Arrange the 1000 results from low to high and select the 10th result as the mutation frequency threshold. SNP sites exceeding the mutation frequency threshold are filtered out;

[0070] S72. Filter out phenotypic outliers in the mutant population: Use R software to draw a box plot, as Figure 2 shown, determine 2 extreme outliers, and retain the phenotypic and genotype data of the remaining 1798 EMS population samples for subsequent analysis;

[0071] S73. Calculate the overall mutation effect score of each sample in each gene based on the genotype annotation information, and obtain the sample genotype mutation effect score matrix. The results are shown in Table 4;

[0072] Table 4 Sample genotype mutation effect scores

[0073]

[0074] Among them, for samples without mutation sites on the gene, the gene mutation effect is recorded as 0 and included in the analysis.

[0075] S74. Based on the generalized linear model, combine the phenotypic data and gene mutation effect score data of each sample, and conduct an association analysis between the genotype mutation effect score and its thousand-grain weight in the R environment, as Figure 3 and Figure 4 shown. Taking the mutation effect score as x and the thousand-grain weight data as y, construct the generalized linear model y = ax + b, perform regression fitting on the model to obtain the fitting linear model coefficients, obtain the significant level of the coefficients through t-test, and obtain the association P-value of the gene;

[0076] S75. Use the box plot drawn by the R software as Figure 2 shown to determine 65 significantly different samples outside the [Q1 - 3IQR, Q3 + 3IQR] interval, where: Q1 and Q3 represent the lower quartile and upper quartile of the thousand-grain weight data of the samples respectively, IQR represents the interquartile range. Group all samples according to the mutation effect score of the gene, with one group having a mutation effect score of 0 and the other group having a non-zero mutation effect score. If the gene is associated with the phenotype, then in the group with a non-zero mutation effect score, the phenotype of the samples should also be significantly different from the overall. Test whether the phenotype of the mutated samples of each gene has a significant difference through chi-square test, and obtain the association P-value of the gene;

[0077] S76. Conduct a permutation test to calculate the significance threshold. After shuffling the correspondence between the sample phenotype and the sample number, perform the calculations in steps S73 and S74, respectively count the minimum values of the P-values, repeat 1000 times, respectively obtain 1000 minimum P-value results for the two methods, sort them from smallest to largest according to the P-values respectively, and select the P-value ranked 50th as the significance threshold. Calculate the -log value of the average of the two p-values, and finally the threshold is determined to be 3.2.

[0078] Take TraesKN3D01HG16320 as an example to illustrate the results generated by the above steps:

[0079] The situation of TraesKN3D01HG16320 is shown in Table 5.

[0080] Table 5 Situation of TraesKN3D01HG16320

[0081]

[0082] The significance threshold for step S74 is 0.000089, and for step S74 is 0.000103.

[0083] The P-value of the association between TraesKN3D01HG16320 and the phenotype calculated using the method of step S73 is 0.000032, ranking 103rd among all genes (a total of 108326 genes).

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

[0085] Therefore, this gene is determined to be a candidate gene related to the thousand-grain weight.

[0086] The list of associated genes in the high-support group based on the new association analysis method is shown in Table 6.

[0087] Table 6 List of associated genes in the high-support group based on the new association analysis method

[0088]

[0089] As can be seen from Table 6, the thousand-grain weight, as a complex trait, is affected by multiple factors. The genes in Table 6 regulate wheat growth and development from multiple perspectives and conform to the developmental characteristics of the population. The method provided by the present invention is applicable to various crop populations mutagenized by EMS, such as rice, wheat, corn, etc. The target trait can be a quantitative trait or a multiple qualitative trait.

[0090] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. A genome-wide association analysis algorithm at the gene level based on the EMS population, characterized in that It includes the following steps: S1. Obtain the sample to be analyzed for exon sequencing; S2. Obtain the phenotypic data of the EMS samples and normal control samples; S3. Align the exon sequencing data of the EMS population with the genome of its species; S4. Use a SNP detection tool to detect SNPs in the obtained alignment results to obtain SNP information in the gene coding region; S5. Use the SnpEff software to annotate the SNPs and predict the mutation effects to obtain the degree of influence of each SNP variation on gene function; S6. Statistically analyze the total mutation effects received in each gene region of each sample to obtain information on the degree of influence at the gene level; S7. Use a linear model and chi-square test for genotype-phenotype association analysis, which specifically includes the following steps: S71. Filter false positive SNPs; S72. Filter phenotypic outliers; S73. Establish a positive statistical model of the mutation effect at the gene level and the phenotype, conduct genotype-phenotype association analysis, and obtain the association P value between the phenotype and the gene; S74. Establish a reverse statistical model of the mutant phenotype and the mutation effect at the gene level, conduct phenotype-genotype association analysis, and obtain the association P value between the gene and the phenotype; S75. Combine the forward and reverse results to obtain the final association score between the phenotype and the gene; S76. According to the principle of statistical test, use permutation testing to obtain a reliable significance threshold; S77. Screen significantly associated genes.

2. The genome-wide association analysis algorithm at the gene level based on the EMS population according to claim 1, wherein In step S4, the SNP information is saved in the VCF file format. The specific method is as follows: Extract the genotype information of the SNP sites of all samples from the VCF file to obtain the genotyping matrix of the SNP sites, and filter out abnormal sites with too high mutation frequencies according to the mutation frequency threshold determined by permutation testing.

3. The gene-level genome-wide association analysis algorithm based on the EMS population according to claim 1, wherein, In step S5, the specific method is as follows: Use the annotation tool SnpEff to annotate the specific gene where each site is located, set the distance threshold K, where K is set as 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 within the gene region, and predict the influence effect of each SNP. The levels of the influence effect are divided into HIGH, MODERATE, LOW, and MODIFIER.

4. The gene-level genome-wide association analysis algorithm based on the EMS population according to claim 1, wherein In step S73, for genotype-phenotype association analysis, it specifically includes the following steps: S731. Assign scores to each site according to the level of the influence effect of the site mutation, where MODIFIER = 1.000001, LOW = 2.00001, MODERATE = 3.0001, HIGH = 4.001; S732. According to the genotype information of the samples and gene position annotation, statistically analyze the total mutation influence effect scores of each gene of each sample; S733. Establish a linear regression statistical model of genotype and phenotype at the gene level, conduct genotype-phenotype association analysis, obtain the association P value between the phenotype and each gene, and select the top 0.5% as significantly associated genes; S734. Establish a test model for the gene mutation situation and the phenotypic situation, conduct an association analysis between gene mutations and phenotypic abnormalities, obtain the association P-values between the phenotype and each gene, and select the top 0.5% as significantly associated genes; S735. Mutually verify the significantly associated genes in steps S3-3 and S3-4, and take the intersection after merging as candidate associated genes; S736. Determine the significance threshold through permutation testing, filter the candidate associated genes with significance less than the threshold, and obtain the final result.

5. The gene-level genome-wide association analysis algorithm based on the EMS population according to claim 1, wherein In step S77, screening significantly associated genes specifically includes the following steps: S771. Set a distance threshold K, where K is set as a positive number less than 5000, and obtain the SNP sites on the gene and within K bp upstream and downstream of the gene; S772. Use simulated data to predict the mutation frequency threshold of SNPs in the population, and filter the SNP sites exceeding the threshold; S773. Conduct a comprehensive calculation based on the association P-values of each gene calculated in steps S73 and S74 to obtain the final association score of the gene; S774. Determine the threshold p0 of the association score, and screen out the genes with an association score greater than p0.

Images