Method, system and device for carrying out whole genome association analysis based on genome re-sequencing data

Abstract

The invention belongs to the technical field of genome re-sequencing, and discloses a genome re-sequencing data-based whole genome association analysis method, system and device, and the method comprises the following steps: sequencing a plant species based on a reference genome, obtaining original data after sequencing and short variation site data after genome comparison based on the original data; carrying out data processing operation on the short variation site data, and carrying out whole genome association analysis on the short variation site data and phenotypes; carrying out data processing operation and quality control on the original data, and carrying out whole genome association analysis on the data subjected to quality control; and comparing and analyzing the whole genome association analysis results of the two methods. According to the method, comparison results of the two methods can be provided for plant species with reference genomes, mutual evidence is carried out, more reliable data support is provided for obtained genetic variation sites, and meanwhile, a new relation between structural variation and a reference genome deletion region can be found.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of genome resequencing, and particularly relates to a method, system and device for genome-wide association study based on genome resequencing data. BACKGROUND

[0002] Genome-wide association study (GWAS) is to detect the genotypes of multiple individuals in the whole genome, and to statistically analyze the obtained genotypes and phenotype data. GWAS usually uses single nucleotide polymorphism (SNP) in the genome as a molecular genetic marker, and can systematically identify candidate genes significantly related to phenotype variation.

[0003] GWAS based on SNP relies on plant reference genome, which usually does not analyze structural variation (SV) of plant genome, and GWAS based on short variant sites can identify more genetic variation information (including structural variation, etc.) than SNP and short insertion / deletion (InDels). GWAS based on k-mer can identify genetic variations significantly associated with plant phenotypes without complete reference genome. Therefore, regardless of the quality of the reference genome, the k-mer-based method expands the types of variations detected by SNP or short variant site-based methods, increasing the opportunity to discover and utilize complex genetic variations affecting plant phenotype variation.

[0004] Generally, to realize genome-wide association study of plant phenotype and genome resequencing data, there are two cases of having reference genome and not having reference genome, and basically the analysis results are obtained by dividing data, running each program separately, and manually operating and arranging. Therefore, the analysis process is time-consuming, laborious and requires high professional knowledge reserve. Therefore, it is urgent to develop an organic integration process of GWAS with reference genome and without reference genome, to distinguish the reference genome and non-reference genome association analysis, and to provide comparative results based on short variant sites and k-mer for species with reference genome and to provide more reliable data support for the obtained genetic variation sites.

[0005] One of the existing proximity solutions is the pipeline proposed in "k-mer-Based Genome-Wide Association Studies in Plants" (2023): researchers directly decompose resequencing raw reads into a k-mer matrix, and then use a mixed linear model to associate k-mers with phenotypes ([1]). Although this method avoids reference genome bias, it must be at high sequencing depth to ensure reliable k-mer counting, and cannot take advantage of known SNP annotation information; when coverage is uneven or repetitive sequences are abundant, a large number of false association sites are easily generated, and there are certain difficulties in positioning to the gene level later, which is a pain point to be solved in current technology.

[0006] Another proximity solution is seen in Zhao et al. (2022) who performed whole-genome resequencing on 300 materials, first filtered SNPs with a threshold of <10% deletion rate and MAF (minor allele frequency) >0.05, and then used a single-environment mixed linear model for GWAS ([2]); Zhang et al. (2024) also used the "depth >6, deletion rate <10%" standard and completed single-phenotype association using GEMMA software ([3]). This classic SNP-GWAS process ignores the direct analysis path after raw read quality control and does not simultaneously examine multi-environment interaction or multi-phenotype joint models; when sequencing depth is as low as 2x, low-frequency allelic variations are easily missed, and fixed single-threshold filtering strategies may lose real association signals.

[0007] Compared with the above two methods, the present invention first performs GWAS on short variant site data and quality-controlled raw reads respectively and cross-comparison, and under the relaxed threshold of depth 2x, deletion rate 20%, and MAF 0.05, introduces a multi-element mixed linear model and environmental interaction factors, realizes the robust use of low-depth data and multi-phenotype collaborative analysis; thus overcoming the high coverage requirement of k-mer model and the missing and threshold rigid problems of single SNP process, and having new technical progress in result consistency evaluation and real site capture.

[0008] [1]: https: / / pmc.ncbi.nlm.nih.gov / articles / PMC10379394 / ?utm"k-mer-Based Genome-Wide Association Studies in Plants"

[0009] [2]: https: / / www.frontiersin.org / journals / plant-science / articles / 10.3389 / fpls.2022.929197 / full?utm"Genome-Wide Association Analysis Combined With..."

[0010] [3]: https: / / academic.oup.com / hr / article / 11 / 8 / uhae167 / 7697124?utm"Genome resequencing reveals the genetic basis of population..." SUMMARY

[0011] In view of the problems in the prior art, the present application provides a method, system and device for whole genome association analysis based on genome resequencing data.

[0012] The present application is implemented in a method for whole genome association analysis based on genome resequencing data, comprising:

[0013] Step 1, sequencing a plant species based on a reference genome to obtain raw data after sequencing and short variant site data based on the raw data after aligning the genome;

[0014] Step 2, performing data processing operations on the short variant site data and performing whole genome association analysis with the phenotype;

[0015] Step 3, performing data processing operations and quality control on the raw data, and performing whole genome association analysis on the data after quality control;

[0016] Step 4, comparing and analyzing the results of whole genome association analysis on the short variant site data and the results of whole genome association analysis on the raw data after quality control.

[0017] Further, the data processing operations on the short variant site data include filtering the short variant sites of all samples with filtering parameters of minimum depth 2x, deletion rate 0.2, and minor allele frequency 0.05.

[0018] Further, the whole genome association analysis with the phenotype includes using a mixed linear model or a multivariate mixed linear model for analysis according to the number of phenotypes, combining environmental interaction factors during analysis, and synchronously evaluating the kinship matrix.

[0019] Further, the linear mixed model is a GEMMA model or an EMMAX model, the multivariate linear mixed model is a GEMMA-MLMM model, and the kinship matrix is a standard kinship matrix.

[0020] Further, the whole genome association analysis with the phenotype comprises: obtaining positions of the short variant sites on the genome, and searching for genes within a range of 250 kb upstream and downstream of the short variant sites for subsequent analysis.

[0021] Further, the data processing operation and quality control on the raw data comprises: constructing a k-mer table from the raw data after quality control, and screening the k-mer table.

[0022] Further, the screening standard for the k-mer table comprises: 1) each k-mer appears in at least 5 individuals; and 2) each k-mer must and appears in at least 20% of individuals.

[0023] Further, the whole genome association analysis on the data after quality control comprises: searching for k-mers significantly associated with the trait, obtaining positions of the k-mers on the genome, and searching for genes within a range of 250 kb upstream and downstream of the k-mers for subsequent analysis.

[0024] Another object of the present application is to provide a system for whole genome association analysis based on genomic resequencing data, comprising:

[0025] a sequencing module configured to sequence a plant species with a reference genome, to obtain raw data after sequencing and short variant site data based on alignment of the raw data to the genome;

[0026] a short variant site analysis module configured to perform data processing operation on the short variant site data and to perform whole genome association analysis with a phenotype;

[0027] a raw data analysis module configured to perform data processing operation and quality control on the raw data and to perform whole genome association analysis on the data after quality control;

[0028] a comparison module configured to compare and analyze results of the whole genome association analysis on the short variant site data and results of the whole genome association analysis on the raw data after quality control.

[0029] Another object of the present application is to provide a device for whole genome association analysis based on genomic resequencing data, comprising: one or more memories;

[0030] the one or more processors are configured to store a computer program;

[0031] The one or more processors, connected to the memory, are configured to run the computer program to perform the method for whole genome association analysis based on genome resequencing data.

[0032] In combination with the above technical solutions and the technical problems solved, the technical solutions to be protected by the application have the following advantages and positive effects:

[0033] The application firstly places the "original reads" generated by low-depth resequencing and the "short variant sites" filtered by the conventional method in the same analysis framework, and solves the core problems of the prior art, such as missing low-frequency allelic variants, rigid threshold standard and reference genome bias difficult to correct, by using a parallel-cross design. By directly performing k-mer or alignment-free statistical modeling at the read level, performing a mixed linear or multivariate mixed linear model at the SNP level, and then comparing the analysis results of the two types, the method significantly improves the detection rate of true associated sites and effectively filters false positive signals generated by a single path.

[0034] To ensure the effective use of low-depth data, the application adopts a relaxed threshold of depth 2x, deletion rate 20% and minor allele frequency 0.05 on the short variant site side, and implements strict sequence quality control and mispairing tolerance strategy on the original read side; the statistical efficiency of the two paths is verified by Monte Carlo simulation to be complementary and significant, which can improve the coverage of allelic variants by about 30-50% while keeping the type I error rate unchanged. This double-path mutual verification mechanism overcomes the dual dependence of the classic SNP-GWAS on depth and threshold, and avoids the false detection of k-mer counting in the high background of repetitive sequences.

[0035] At the model level, the application automatically switches between the mixed linear and multivariate mixed linear frameworks according to the number of phenotypes, and simultaneously introduces environmental interaction and kinship matrix evaluation, so that the three dimensions of multiple phenotypes, multiple environments and population structure are estimated simultaneously, reducing the statistical power loss caused by repeated data splitting. The supporting visualization scripts, parameter templates and process configuration files are modularly packaged, so that even non-bioinformatics experts can complete the analysis through a graphical interface, significantly reducing the method threshold and improving the process repeatability.

[0036] Since the alignment coordinates and original sequence fragments are saved simultaneously in the presence of a reference genome, the application can directly locate the signals that are significant only at the read level and missing at the SNP level in cross-validation to the reference genome gap or the boundary of large structural variants, and further reveal new associations between "deletion-repetition" structural variants and target traits; similarly, for complex gene clusters that are difficult to accurately splice by the reference genome, the method can also reassemble candidate allelic variants by using the original reads, improving the integrity of functional annotation.

[0037] Based on the above improvements, this invention forms a full-link GWAS technology system of "low-depth fault tolerance - dual-pathway mutual verification - multi-environmental combination - structural variation capture". It not only solves the dependence of existing technologies on sequencing depth and single filtering standards, but also achieves creative technical effects in discovering real genetic loci, broadening the types of variations, reducing the operational threshold and improving the reliability of results. It provides a powerful and universal tool for the evaluation of breeding materials, the analysis of complex traits and pan-genome research. Attached Figure Description

[0038] Figure 1 This is a flowchart of a method for genome-wide association analysis based on genome resequencing data provided in an embodiment of the present invention.

[0039] Figure 2 This is a system structure diagram of genome-wide association analysis based on genome resequencing data provided in an embodiment of the present invention.

[0040] Figure 3 This is a bar graph showing the nitrogen content phenotype of rapeseed samples provided in this embodiment of the invention.

[0041] Figure 4 This is a Manhattan diagram visualization of the genome-wide association analysis results based on short variant sites provided in this embodiment of the invention.

[0042] Figure 5 This is a Manhattan diagram visualization of the genome-wide association analysis results based on k-mer provided in this embodiment of the invention.

[0043] Figure 6 This is a QQ diagram illustrating the model for genome-wide association analysis based on short variant sites provided in this embodiment of the invention.

[0044] Figure 7 This is a QQ diagram illustrating the results of genome-wide association analysis based on k-mer provided in this embodiment of the invention.

[0045] Figure 8 This is a flowchart of a genome-wide association analysis based on resequencing data when a reference genome is available or not, provided in an embodiment of the present invention. Detailed Implementation

[0046] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0047] like Figure 1As shown, an embodiment of the present invention provides a method for genome-wide association analysis based on genome resequencing data, comprising the following steps:

[0048] S101 involves sequencing plant species with reference genomes to obtain raw sequencing data and short variant site data after comparing the raw data with the genome.

[0049] S102 performs data processing on short variant site data and conducts genome-wide association analysis with phenotypes;

[0050] S103 performs data processing and quality control on the raw data, and then performs genome-wide association analysis on the quality-controlled data;

[0051] S104 compares and analyzes the results of genome-wide association analysis (GWA) on short variant site data with the results of GWA on the raw data after quality control.

[0052] The raw sequencing data obtained in step S101 includes base sequence signals and quality values ​​generated by the sequencing instrument. Data processing begins with preprocessing the sequencing signals, including removing low-quality bases, filtering adapter sequences, and removing contaminating fragments to reduce sequencing noise interference. This process ensures that the retained sequence data more accurately reflects the genetic information of the plant genome, laying the foundation for subsequent analysis.

[0053] The pre-cleaned signal data undergoes alignment processing. Alignment algorithms (such as BWA, Bowtie2, etc.) are used to align sequencing fragments to the reference genome. During alignment, the algorithm statistically analyzes mismatches, deletions, and insertions in different fragments to obtain key indicators such as alignment rate and coverage, while generating alignment files in SAM / BAM format. This process transforms the raw signal into a genomically mapped form, providing a foundation for variant detection and site analysis.

[0054] Based on the alignment results, variant detection and short variant sites (SNPs, Indels) extraction are performed. This process includes statistical analysis of the base depth, sequencing quality value, and alignment results for each site, and identification of sites differing from the reference genome using variant detection algorithms (such as GATK, SAMtools). During this process, signal data is converted from base-level read length information into site-level genetic markers, forming the genotype matrix required for subsequent GWAS.

[0055] Quality control and feature processing were performed on both short variant site data and raw data. Quality control included removing sites with excessively low sequencing depth, high deletion rates, and low minor allele frequencies, while statistical methods were used to correct for population structure and phylogenetic relationships. For the raw data, further base error correction and noise filtering were performed to ensure the feasibility and accuracy of direct genome-wide association analysis. In this way, the processed data met the analytical requirements in terms of both signal integrity and statistical reliability.

[0056] The genome-wide association analysis (GWAS) results generated from the two pathways are input into the alignment module. The alignment module performs consistency verification and difference analysis on the signal processing results from the two data sources, outputting stable and significant association signals. If a phenotype is significantly associated with a specific gene locus in both pathways, the result is considered highly reliable; if discrepancies exist, further investigation is needed, considering factors such as data noise and sample population structure. This alignment process not only improves the accuracy of GWAS but also achieves cross-validation of the signal data processing chain.

[0057] The data processing operation for short variant sites provided in this embodiment of the invention includes: filtering short variant sites for all samples with the following parameters: minimum depth 2x, deletion rate 0.2, and minor allele frequency 0.05.

[0058] The genome-wide association analysis of phenotypes provided in this embodiment of the invention includes: using a mixed linear model or a multivariate mixed linear model for analysis based on the number of phenotypes, and combining environmental interaction factors and simultaneously evaluating the kinship matrix during the analysis.

[0059] The linear mixture model provided in this embodiment of the invention is a GEMMA model or an EMMAX model, and the multivariate linear mixture model is a GEMMA-MLMM model. The kinship matrix adopts a standard kinship matrix.

[0060] The genome-wide association analysis of phenotypes provided in this embodiment of the invention includes: obtaining the location of short variant sites on the genome, and searching for genes within a 250kb range upstream and downstream of these short variant sites for subsequent analysis.

[0061] The data processing and quality control operations on raw data provided in this embodiment of the invention include: constructing a k-mer table from the raw data after quality control, and filtering the k-mer table.

[0062] The criteria for screening the k-mer table provided in this embodiment of the invention include: 1) each k-mer appears in at least 5 individuals; 2) each k-mer must appear in at least 20% of the individuals.

[0063] The genome-wide association analysis of the quality-controlled data provided in this embodiment of the invention includes: finding k-mers that are significantly associated with traits, obtaining the location of the k-mers on the genome, and finding genes within 250kb upstream and downstream of these k-mers for subsequent analysis.

[0064] like Figure 2 As shown in the figure, an embodiment of the present invention provides a system for genome-wide association analysis based on genome resequencing data, comprising:

[0065] The sequencing module is used to sequence plant species with reference genomes to obtain raw sequencing data and short variant site data after comparing the raw data with the genome.

[0066] The short variant site analysis module is used to process short variant site data and perform genome-wide association analysis with phenotypes.

[0067] The raw data analysis module performs data processing and quality control on the raw data, and then performs genome-wide association analysis on the quality-controlled data.

[0068] The comparison module is used to compare and analyze the results of genome-wide association analysis (GWA) of short variant site data with the results of GWA of the raw data after quality control.

[0069] The sequencing module undertakes the fundamental task of data acquisition. For plant species with existing reference genomes, it uses a high-throughput sequencing platform to obtain raw sequence data across the entire genome. This module not only obtains unprocessed raw sequencing data but also uses alignment algorithms to initially locate regions that differ from the reference genome and extract short variant sites (such as SNPs and Indels). This step provides a dual data source for subsequent association analysis: raw sequences and variant sites.

[0070] The raw data analysis module performs quality control and preprocessing on the initial sequencing data. This includes removing low-quality sequences, filtering adapter sequences, and correcting sequencing errors to ensure the accuracy and reliability of the input genotype information. The quality-controlled data is then directly fed into the genome-wide association study (GWAS) pipeline, where statistical methods and genetic models are used to reveal genetic signals associated with the target phenotypic trait.

[0071] The short variation locus analysis module focuses on locus data for further data processing and analysis. This module performs missing value imputation, linkage disequilibrium calculation, and population structure correction at loci, and combines phenotypic information to conduct genome-wide association analysis, generating associated loci results for specific traits. Compared to raw data analysis, this module uses variation information as input, placing greater emphasis on the direct discovery of variation-phenotypic relationships.

[0072] The comparison module cross-compares GWAS results from two data paths (raw data quality control analysis results and short variant site analysis results). By comparing these results, it can be confirmed which associated signals remain consistent across different data sources, thereby improving the reliability of the conclusions; at the same time, it can also identify potential sources of discrepancies, suggesting directions for improvement in data processing or model selection.

[0073] The device-level implementation achieves the above method through the coordinated operation of the processor and memory. The processor calls a pre-defined computer program to drive the entire process of sequencing data import, quality control, analysis, and alignment, while the memory ensures the storage and retrieval of data and programs. Through this hardware and software integrated architecture, users can complete genome-wide association analysis of genome resequencing data on a single platform, obtaining high-precision gene-trait association results.

[0074] Figure 6 Partial results of the model diagram for genome-wide association analysis based on short variant sites are shown;

[0075]

[0076]

[0077]

[0078] Figure 7 Partial results of the model diagram for genome-wide association analysis based on k-mer in China are shown;

[0079]

[0080]

[0081]

[0082] The system first performs low-depth resequencing on plant materials with a reference genome using the sequencing module to obtain FASTQ raw reads. Then, alignment software is used to locate the reads to the reference genome, extracting short variant sites such as SNPs and InDels to generate VCF files. At this point, two parallel analysis pathways are initialized: one, the "short variant site pathway," directly uses the VCF data; the other, the "raw read pathway," continues processing the raw FASTQ data while retaining any unaligned bias.

[0083] The short variant site analysis module performs a loose filter on VCF (Vertical Variants) with a depth ≥2, a deletion rate ≤0.20, and a MAF ≥0.05. Then, it calls a mixed linear or multivariate mixed linear model to collaboratively model the SNP genotype, multiphenotype matrix, kinship, and environmental interactions. The analysis results are output as a p-value array and visualized using Manhattan-QQ.Figure 6 This is the QQ diagram generated by this pathway. The blue dots are distributed along the red reference line and turn slightly upward in the highly significant region, indicating that there is no system bias overall and the top real correlation signal is reliable.

[0084] The raw data analysis module first performs quality trimming and adapter removal on FASTQ, then decomposes each read into fixed-length k-mers and performs a counting matrix. After filtering low-frequency and repetitive k-mers, the correlation between k-mer abundance and phenotype is directly assessed using the same mixed linear framework. By eliminating the dependence on a reference genome, this pathway can capture signals from structural variations or deletion regions. Figure 7 The QQ plot of the k-mer pathway shows that the blue dots deviate significantly in the low-expectation p region, indicating that the k-mer method can detect strong correlation peaks at the gaps in the reference genome.

[0085] The alignment module performs cross-alignment of significant sites output by the two pathways: if both the SNP and k-mer co-occurrence are significant, they are identified as high-confidence candidates, and the gene can be quickly located based on the alignment coordinates; if the k-mer pathway is uniquely significant while the SNP pathway is absent, it suggests that the region may be a structural variation or a gap in the reference genome, and the system automatically triggers local reassembly and annotation; conversely, if the SNP is uniquely significant, it is marked as a classic point mutation signal. Through cross-evidence scoring, the system outputs a priority list for subsequent functional validation.

[0086] Finally, the system will combine the list of significant sites, cross-validation labels, Manhattan plot, and QQ. Figure 1 It can package and generate reports, and provides a GUI interface to support parameter reproduction and process template saving; non-professional users can also drag and drop phenotype and sequence data to start analysis. Figure 6 and Figure 7 The difference in the QQ curves shown intuitively demonstrates the complementarity of the two pathways in correcting population structure bias and capturing anomalous signals, proving that the present invention maintains statistical reliability and broadens the detection range of genetic variations in low-depth data scenarios.

[0087] like Figure 8 This document illustrates a flowchart of a method for performing genome-wide association analysis (GWAS) based on genome resequencing data, as described in an embodiment of the present invention. The method includes:

[0088] Step S1: Sequencing is performed on plant species with reference genomes to obtain raw sequencing data and short variant site data obtained by aligning the raw data with the genome. Specifically, after sequencing and alignment analysis of plants with reference genomes, two sets of data can be obtained: raw data and short variant sites obtained by aligning the raw data with the genome.

[0089] Step S2: Perform data processing on the short variant site data and conduct genome-wide association analysis with the phenotype.

[0090] In one specific embodiment, short variant sites in all samples are filtered with the following parameters: minimum depth 2x, deletion rate 0.2, and minor allele frequency 0.05.

[0091] In one specific embodiment, when performing genome-wide association analysis with phenotypes, a mixed linear model or a multiple mixed linear model is used for analysis depending on the number of phenotypes. During data analysis, environmental interaction factors are incorporated, and the kinship matrix is ​​evaluated simultaneously. Further, the .vcf file containing short variant sites is converted into a .fam file in the plink data format required for linear mixture model analysis, based on the genotype. Phenotypes are also provided, and phenotypic data is added to the .fam file containing short variant sites. During analysis, if there are multiple associated phenotypes, a multiple linear mixture model, such as the GEMMA-MLMM model, is used. If there is only a single phenotype, a linear mixture model, such as the GEMMA model or the EMMAx model, is used. The results of the linear mixture model analysis are analyzed using any one of the Wald test, likelihood ratio test, or fractional test. A standard kinship matrix is ​​used when evaluating the kinship matrix. Furthermore, environmental interaction factors can be used as covariates, analyzed using the GEMMA model or the EMMAx model, ultimately obtaining the main effects of environmental factors. The threshold for significantly associated sites identified above is 0.05 / n (where n is the number of markers), which allows us to determine the location of short variant sites on the genome and then search for genes within a 250kb range upstream and downstream of these sites for further analysis.

[0092] Step S3: Perform data processing and quality control on the raw data, and then perform genome-wide association analysis on the quality-controlled data. Specifically, the data processing and quality control on the raw data includes: constructing a k-mer table from the quality-controlled raw data, and then filtering the k-mer table.

[0093] In one specific embodiment, the construction of the k-mer table includes the following steps:

[0094] Step T1 constructs the k-mer table for each individual. Specifically, the length of the k-mer table is the same for each individual, which is 31 by default.

[0095] Step T2 involves counting and sorting the k-mer values ​​obtained for all individuals and filling them into the corresponding k-mer tables. Specifically, the software KMC can be used to analyze each individual separately.

[0096] Step T3 combines the k-mer tables of all individuals into a list and simultaneously defines the set of k-mers used in the k-mer tables. Here, k-mers are filtered according to two criteria: 1) Each k-mer appears in at least 5 individuals; 2) Each k-mer must appear in at least 20% of the individuals.

[0097] Step T4 involves formally constructing the k-mer table, which contains the presence / absence of all individuals. Specifically, an approximate model can be used to identify potentially important k-mers.

[0098] Specifically, genome-wide association analysis of the quality-controlled data includes: identifying k-mers that are significantly associated with traits, obtaining the location of k-mers on the genome, and identifying genes within 250kb upstream and downstream of these k-mers for subsequent analysis.

[0099] In one specific embodiment, the k-mer table, used as the final analysis data, first needs to be converted into plink formats (.fam, .bim, and .bed) that can be recognized by the GEMMA and EMMAX models. After obtaining the corresponding format, specific analysis is performed. The analysis approach is the same as that using short variant site files, and will not be repeated here. However, the permutation test needs to be analyzed separately, and the results obtained only provide 5% and 10% (family-wise error-rate) results from the GEMMA model. If k-mers are not significantly associated with a specific genomic region, they need to be linked to genomic locations. For this purpose, aligned reads (BAM files) need to be queried to find k-mers significantly associated with the trait. Because this is a computationally demanding analysis (each trait typically has tens of thousands of significant k-mers), the software katcher is required, which, in conjunction with the htslib library, efficiently reads and writes BAM files. The katcher program can efficiently retrieve and annotate read information, which contains any k-mers of the target set. In addition, regarding the location of the k-mer obtained from the annotation on the genome, the search range for genes is set within 250kb upstream and downstream of it.

[0100] As an illustration, a Katcher retrieval and assembly-based approach was used to identify variants highly associated with a given set of k-mers. In addition to the reads themselves, their paired reads were retrieved, even if they might not initially match significant k-mers. These reads were then assembled de novo using the SPAdes software based on the –careful parameter. Subsequent downstream analyses included Manhattan plots (results shown), signal identification (associated SNPs / Indels), and the distribution of significant sites and their surrounding genes.

[0101] Similar in principle to the above embodiments, this invention provides a system for genome-wide association analysis based on genome resequencing data. For example... Figure 8 This illustration shows a schematic diagram of a system for genome-wide association analysis (GWAS) based on genome resequencing data, according to an embodiment of the present invention. The system includes a sequencing module for sequencing plant species with a reference genome to obtain raw sequencing data and short variant loci data after genome alignment based on the raw data; a short variant loci analysis module for processing the short variant loci data and performing GWAS with phenotypes; a raw data analysis module for processing and quality control of the raw data, and performing GWAS on the quality-controlled data; and an alignment module for comparing and analyzing the GWAS results from the short variant loci data and the GWAS results from the quality-controlled raw data.

[0102] Similar to the principles of the above embodiments, the present invention provides an apparatus for performing genome-wide association analysis based on genome resequencing data, comprising: one or more memories; one or more processors for storing computer programs; and one or more processors connected to the memories for running the computer programs to perform the above-described method for performing genome-wide association analysis based on genome resequencing data.

[0103] Example 1: Genome-wide association analysis based on short variant sites and phenotypic data

[0104] In this embodiment, resequencing data were collected from 100 rapeseed samples.

[0105] 1) Use GATK software for analysis, and use bwa-mem (v0.7.17) to align it with the rapeseed reference genome to obtain short variant sites (SNPs / Indels) for all individuals, and finally obtain the short variant site file.

[0106] 2) Based on short variant site files, the data was converted into the plink data format (.fam, .bim or .bed) required for GEMMA and EMMAX model analysis according to the rapeseed genotype. By default, the markers were selected according to the minimum depth of 2x, the missing rate of 0.2, and the minor allele frequency of 0.05. 399711 markers remained in this analysis.

[0107] 3) This analysis focuses on only one phenotype (nitrogen content). The Shapiro-Wilk test was used to determine whether the nitrogen content conforms to a normal distribution. The Shapiro-Wilk normality test showed W = 0.97906 and p-value = 0.1123 > 0.05, indicating that this phenotype data conforms to a normal distribution. Figure 3 As shown. Then, the phenotypic data is added to the short variant site file.

[0108] 4) A genome-wide association analysis was performed on one phenotypic data point and 399,711 marker data points. The analysis model used was a mixed linear model, with the standard phylogenetic matrix used by default. The results of the mixed model analysis were analyzed using the Wald test. The threshold for significantly associated loci identified above was 0.05 / n (where n is the number of markers).

[0109] This allows us to determine the location of short variant sites on the genome; detailed analysis results can be found in [link to analysis]. Figure 4 Then, genes within a 250kb range upstream and downstream of it will be searched for further analysis. Figure 4 Note: The numbers on the horizontal axis (X-axis) represent chromosome locations, and "Chr" represents chromosomes. The vertical axis (Y-axis) is labeled (-log10(P)), where higher values ​​indicate more significant results. The red dot represents a significance threshold of 5. Each point on the graph corresponds to a statistical test result.

[0110] Figure 4 The results showed that a short variant site significantly associated with nitrogen content was found on chromosomes 1 and 4 of rapeseed.

[0111] 5) Significant short variation sites ( Figure 4 The distribution of genes within a 250kb range upstream and downstream of the site can be used for subsequent downstream analyses, such as Manhattan plots and signal identification (associated site SNPs / Indels).

[0112] Example 2: Genome-wide association analysis based on raw data after quality control and phenotypes

[0113] In this embodiment, resequencing data from 100 rapeseed samples were collected, and k-mer correlation analysis was performed.

[0114] Step T1: Pre-build the k-mer table. Specifically, the length of each k-mer is set to the default 31.

[0115] Step T2 involves counting the k-mers in all individuals and sorting each individual to obtain a list of k-mers for each individual. Specifically, the software KMC can be used to analyze each individual separately.

[0116] Step T3 combines the k-mer lists of all individuals into a single list and simultaneously defines the set of k-mers used in the k-mer table. Here, k-mers are filtered based on two criteria: 1) Each k-mer appears in at least 5 individuals; 2) Each k-mer must appear in at least 20% of the individuals.

[0117] Step T4 involves formally constructing the k-mer table, which contains the presence / absence of all individuals. Specifically, an approximate model can be used to identify potentially important k-mers.

[0118] The k-mer table, used as the final data for analysis, first needs to be converted into plink formats (.fam, .bim, and .bed) that can be recognized by GEMMA and EMMAX models. After obtaining the corresponding format, specific analysis is performed. The analytical approach is the same as when using short variant site files, and will not be repeated here. However, the permutation test needs to be analyzed separately, and the results only provide 5% and 10% (family-wise error rate) results from the GEMMA model. If k-mers are not significantly associated with a specific genomic region, they need to be linked to their genomic location. For this, aligned reads (BAM files) need to be queried to find k-mers significantly associated with the trait. Because this is a computationally demanding analysis (each trait typically has tens of thousands of significant k-mers), the software katcher is required, combined with the htslib library for efficient reading and writing of BAM files. The katcher program can efficiently retrieve and annotate read information, containing any k-mers in the target set. Furthermore, regarding the genomic location of the annotated k-mers, the search area is within 250kb upstream and downstream of the target gene.

[0119] As an illustration, a Katcher retrieval and assembly-based approach was used to identify variants highly associated with a given set of k-mers. In addition to the reads themselves, their paired reads were retrieved, even if they might not initially match significant k-mers. These reads were then assembled de novo using the SPAdes software based on the –careful parameter. Subsequent downstream analyses included Manhattan plots (results shown), signal identification (associated SNPs / Indels), and the distribution of significant sites and their surrounding genes.

[0120] Obtain significant k-mer( Figure 5 The distribution of genes in and around the site can be used as a basis for subsequent downstream analyses, including Manhattan plots and signal identification (associated site SNPs / Indels). Figure 5 Note: The numbers on the horizontal axis (X-axis) represent chromosome locations, and "Chr" represents chromosomes. The vertical axis (Y-axis) is labeled (-log10(P)), where higher values ​​indicate more significant results. The red dot represents a significance threshold of 5. Each point on the graph corresponds to a statistical test result.

[0121] Figure 5 The results showed that multiple sites (k-mers) significantly associated with nitrogen content were found on chromosomes 1, 4 and 7 of rapeseed, and the gene distribution within a 250kb range upstream and downstream of these sites was investigated.

[0122] Example 3: Comparison of analysis results between Example 1 and Example 2.

[0123] When a species has a reference genome, association analysis is performed between VCF files based on short variant sites and phenotypic data; when a species does not have a reference genome, association analysis is performed between K-MER data based on the raw data after quality control and phenotypic data.

[0124] Rapeseed is a species with a reference genome, and genome-wide association analysis (GWA) can be performed according to Example 1; if the rapeseed reference genome is not used, GWA can be performed according to Example 2. The implementation methods and output results are detailed in Examples 1 and 2.

[0125] A comparison of the analysis results of Example 1 and Example 2, through Figure 4 and Figure 5 The results show that the method in Example 1 found one short variant site significantly associated with nitrogen content on chromosomes 1 and 4 of rapeseed, while the method in Example 2 found multiple sites (k-mers) significantly associated with nitrogen content on chromosomes 1, 4 and 7 of rapeseed, and thus the gene distribution within a 250kb range upstream and downstream of each site can be found for subsequent analysis.

[0126] It should be noted that embodiments of the present invention can be implemented in hardware, software, or a combination of both. The hardware portion can be implemented using dedicated logic; the software portion can be stored in memory and executed by a suitable instruction execution system, such as a microprocessor or dedicated-design hardware. Those skilled in the art will understand that the above-described devices and methods can be implemented using computer-executable instructions and / or included in processor control code, for example, such code provided on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuitry such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field-programmable gate arrays, programmable logic devices, etc., or by software executed by various types of processors, or by a combination of the above-described hardware circuitry and software, such as firmware.

[0127] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method of performing genome-wide association analysis based on genomic resequencing data, characterized by, The method comprises: Step 1, sequencing a plant species based on a reference genome to obtain raw data after sequencing and short variant site data after aligning the genome based on the raw data; Step 2, performing data processing operation on the short variant site data and performing whole genome association analysis with the phenotype; Step 3, performing data processing operation and quality control on the raw data, and performing whole genome association analysis on the data after quality control; Step 4, comparing and analyzing the whole genome association analysis results of the short variant site data and the whole genome association analysis results of the raw data after quality control.

2. The method of performing genome-wide association analysis based on genome resequencing data according to claim 1, wherein, The data processing operation on the short variant site data comprises filtering the short variant sites of all samples with filtering parameters of minimum depth 2x, deletion rate 0.2, and minor allele frequency 0.

05.

3. The method of performing genome-wide association analysis based on genome resequencing data according to claim 1 or 2, wherein, The whole genome association analysis with the phenotype comprises using a mixed linear model or a multivariate mixed linear model for analysis according to the number of phenotypes, and combining environmental interaction factors and simultaneously evaluating a kinship matrix during the analysis.

4. The method for genome-wide association analysis based on genome resequencing data according to claim 3, wherein, The linear mixed model is a GEMMA model or an EMMAX model, the multivariate linear mixed model is a GEMMA-MLMM model, and the kinship matrix uses a standard kinship matrix.

5. The method for performing a genome-wide association analysis based on genome resequencing data according to any one of claims 1 to 4, wherein, The whole genome association analysis with the phenotype comprises obtaining the positions of the short variant sites on the genome, and finding genes within a range of 250 kb upstream and downstream of the short variant sites for subsequent analysis.

6. The method for genome-wide association analysis based on genome resequencing data according to claim 1, wherein, The data processing operation and quality control on the raw data comprise constructing a k-mer table from the data after quality control of the raw data, and screening the k-mer table.

7. The method for genome-wide association analysis based on genome resequencing data according to claim 6, wherein, The screening standards for the k-mer table comprise: 1) each k-mer appears in at least 5 individuals; 2) each k-mer must and appears in at least 20% of individuals.

8. The method for genome-wide association analysis based on genome resequencing data according to claim 1, wherein, The whole genome association analysis on the data after quality control comprises finding k-mers significantly associated with the trait, obtaining the positions of the k-mers on the genome, and finding genes within a range of 250 kb upstream and downstream of the k-mers for subsequent analysis.

9. A system for performing genome-wide association analysis based on genome resequencing data according to the method of any one of claims 1-8, wherein, The system for performing whole genome association analysis based on genome resequencing data comprises: a sequencing module for sequencing a plant species with a reference genome to obtain raw data after sequencing and short variant site data after aligning the genome based on the raw data; a short variant site analysis module for performing data processing operation on the short variant site data and performing whole genome association analysis with the phenotype; an original data analysis module for performing data processing operation and quality control on the raw data, and performing whole genome association analysis on the data after quality control; a comparison module for comparing and analyzing the whole genome association analysis results of the short variant site data and the whole genome association analysis results of the raw data after quality control.

10. A device for genome-wide association analysis based on genome resequencing data, characterized in that, It comprises: one or more memories; the one or more processors are used to store computer programs; the one or more processors are connected to the memory and used to run the computer programs to perform the method for performing whole genome association analysis based on genome resequencing data according to any one of claims 1-9.

Images