Genome hierarchical feature screening method and device and electronic equipment

Abstract

The invention provides a genomic hierarchical feature screening method and device and electronic equipment, and belongs to the technical field of computers, medicine and breeding, and the method comprises the following steps: carrying out chromosome grouping on genotype information on each SNP site, carrying out missing genotype filling on each chromosome after grouping, and carrying out genotype detection; carrying out chromosome merging on all the filled chromosomes to obtain a primary characteristic gene pool; calculating the deviation between the joint allele frequency between the two SNP sites on the whole chromosome in the primary characteristic gene pool and the expected frequency, calculating a linkage imbalance coefficient based on the deviation, and reserving the SNP sites with the linkage imbalance coefficient smaller than a first set value to obtain a secondary characteristic gene pool; and calculating a Pearson's correlation coefficient between each pair of SNP sites in the secondary characteristic gene pool, and reserving the SNP sites with the Pearson's correlation coefficients smaller than a second set value to obtain a tertiary characteristic gene pool. According to the method, the prediction accuracy of the breeding value is improved through three-stage feature screening.

Description

TECHNICAL FIELD

[0001] The present application relates to the fields of computer, medicine and breeding technology, and in particular, to a genomic hierarchical feature screening method and device. BACKGROUND

[0002] Genomic Selection (GS) is a modern genetic evaluation method that uses high-density single nucleotide polymorphism (SNP) marker information within the whole genome range, combined with phenotype data to construct a prediction model, to achieve accurate estimation of individual genomic breeding value (GEBV). This technology significantly improves the accuracy of genetic evaluation and breeding efficiency by integrating whole genome markers and phenotype association, and has become a core means in animal and plant breeding. In this process, feature gene screening as a key preprocessing link aims to identify a subset of representative and high genetic contribution sites from a large number of SNPs, laying a foundation for building a robust prediction model.

[0003] However, the existing feature screening strategy still has obvious defects, which seriously restricts the further improvement of the prediction performance of GS. Current methods mostly rely on a single technical route, each with its limitations: Principal Component Analysis (PCA) may filter out some sites with subtle genetic effects in the dimensionality reduction process; Linkage Disequilibrium (LD) screening is easily affected by population structure, leading to the loss of key markers; Genome-Wide Association Study (GWAS) has insufficient ability to identify micro-effect sites due to the severe multiple testing problem; and Epistasis (Epi) screening faces the dilemma of high computational complexity and low statistical efficiency. In addition, existing prediction models usually use equal weight processing for all markers, ignoring the heterogeneity of different site effect values, further causing information loss and signal dilution, limiting the prediction accuracy.

[0004] In summary, the main problem of the existing technology is that the genomic feature screening strategy is too single and lacks an adaptive weight mechanism, making it difficult to balance coverage and specificity, resulting in the loss of key genetic information. Therefore, developing a feature screening method that can effectively integrate multi-dimensional information, reduce redundancy and retain high-value signals has become an urgent need to improve the accuracy of genomic selection breeding. SUMMARY

[0005] Therefore, it is necessary to provide a genomic hierarchical feature screening method, device and electronic equipment to solve the technical problem that the existing genomic feature screening strategy is too single and lacks an adaptive weight mechanism, resulting in the loss of key genetic information.

[0006] To solve the above technical problems, in a first aspect, the present application provides a genomic hierarchical feature screening method, comprising: obtaining original genotype data, the original genotype data comprising genotype information at each SNP site on a genome and sample information, the genotype information at each SNP site comprising a chromosome number and a SNP site number; grouping the genotype information at each SNP site by chromosome, filling in missing genotypes for each chromosome after grouping, and merging all chromosomes after filling in to obtain a first characteristic gene pool; calculating a deviation between a combined allele frequency and an expected frequency between two SNP sites on an entire chromosome in the first characteristic gene pool, calculating a linkage disequilibrium coefficient based on the deviation, and retaining SNP sites with linkage disequilibrium coefficients less than a first set value to obtain a second characteristic gene pool; calculating a Pearson correlation coefficient between each pair of SNP sites in the second characteristic gene pool, and retaining SNP sites with Pearson correlation coefficients less than a second set value to obtain a third characteristic gene pool.

[0007] In a possible implementation, before grouping the genotype information at each SNP site by chromosome, the method further comprises: calculating a minimum allele frequency of each SNP site, and excluding SNP sites with minimum allele frequencies less than a third set value; performing a balance test on the each SNP site, and excluding SNP sites with significance probability values less than a fourth set value in the balance test; calculating a missing rate of the each SNP site, and excluding SNP sites with missing rates greater than a fifth set value; retaining all SNP sites that are not excluded to obtain filtered genotype data.

[0008] In a possible implementation, the minimum allele frequency of each SNP site is obtained by the following calculation formula:

[0009]

[0010] In the formula, is the minimum allele frequency at the SNP site, , respectively represent the frequencies of A, B alleles, and satisfy , N represents a total number of individuals, and AA, BB, and AB represent three possible genotypes at a SNP site, and a number of AA homozygote genotypes, a number of BB homozygote genotypes and a number of AB heterozygote genotypes, respectively; In a possible implementation, the balance test on the each SNP site comprises: calculating expected frequencies of the three genotypes on the each SNP site, an expected frequency of AA homozygote genotype being , an expected frequency of BB homozygote genotype being , and an expected frequency of AB heterozygote genotype being ; calculating expected numbers of the three genotypes on the each SNP site, the expected numbers being the expected frequencies multiplied by the total number of individuals; calculating a statistic of a chi-square test of the actual observed numbers of the three genotypes and the expected numbers χ 2 , obtained by the following formula:

[0011] In the formula, , and are the expected numbers of the three genotypes; calculating degrees of freedom of the chi-square test, the degrees of freedom being a difference between a number of genotypes and a number of alleles; inquiring a chi-square distribution critical value table by the statistic of the chi-square test and the degrees of freedom to obtain a significance probability value of the chi-square test.

[0012] In a possible implementation, the sample information comprises a phenotype value; the genotype information on the each SNP site is grouped by chromosomes, each chromosome after grouping is filled with a missing genotype, and all chromosomes after filling are merged to obtain a primary characteristic gene pool, comprising: grouping each SNP site in the screened genotype data according to chromosome numbers and saving as a separate file for each chromosome; estimating possible genotypes of missing SNP sites in the each chromosome through linkage blocks in the screened genotype data, and filling genotypes of the missing SNP sites in the chromosome; merging the filled chromosomes to obtain a complete genotype matrix, and performing whole genome correlation analysis according to the phenotype value, screening SNP sites with a significance probability value less than a sixth set value in the whole genome correlation analysis to obtain the primary characteristic gene pool; The linkage block refers to SNP sites closely arranged in a chromosome.

[0013] In a possible implementation, the calculation of the deviation between the joint allele frequency and the expected frequency between two SNP sites on the whole chromosome in the primary feature gene pool comprises: calculating the deviation between the joint allele frequency and the expected frequency , and the linkage disequilibrium coefficient is calculated based on the deviation, comprising:

[0014] In the formula, is the frequency of the haploid genotype A and B, is the frequency of the allele A, is the frequency of the allele B; the deviation between the joint allele frequency and the expected frequency , the linkage disequilibrium coefficient between two SNP sites on each chromosome is calculated , and the linkage disequilibrium coefficient is calculated by the following formula:

[0015] In the formula, is the frequency of the non-A allele of the first SNP site, is the frequency of the non-B allele of the second SNP site.

[0016] In a possible implementation, the calculation of the Pearson correlation coefficient between each pair of SNP sites in the secondary feature gene pool comprises: calculating the Pearson correlation coefficient between each pair of SNP sites, and the Pearson correlation coefficient is calculated by the following formula:

[0017] In the formula, and are the gene encoding values of the first individual at the corresponding SNP sites, i and are the average gene encoding values of the two SNP sites, respectively, is the total number of individuals. N

[0018] In another aspect, the present application also provides a gene screening effect evaluation method, comprising: determining a tertiary feature gene pool based on the genomic hierarchical feature screening method according to any one of claims 1-7; predicting the unknown phenotype values in all SNP sites according to the tertiary feature gene pool and the genotype-phenotype value mapping relationship, to obtain predicted all phenotype values; the genotype-phenotype value mapping relationship is used to represent the correspondence between the tertiary feature gene pool and the phenotype values corresponding to the SNP sites in the tertiary feature gene pool.​ a correlation coefficient between the predicted overall phenotype value and the actual overall phenotype value is calculated, and a gene screening effect is evaluated by the correlation coefficient.

[0019] In a second aspect, the present application further provides a genomic hierarchical feature screening device, comprising: an original genotype data acquisition unit configured to acquire original genotype data, wherein the original genotype data comprises genotype information at each SNP site on a genome and sample information, and the genotype information at each SNP site comprises a chromosome number and a SNP site number; a primary feature gene screening unit configured to group the genotype information at each SNP site by chromosomes, fill in missing genotypes for each chromosome after grouping, and combine all chromosomes after filling in to obtain a primary feature gene pool; a secondary feature gene screening unit configured to calculate a deviation between a combined allele frequency and an expected frequency between two SNP sites on an entire chromosome in the primary feature gene pool, calculate a linkage disequilibrium coefficient based on the deviation, and retain SNP sites with the linkage disequilibrium coefficient less than a first set value to obtain a secondary feature gene pool; a tertiary feature gene screening unit configured to calculate a Pearson correlation coefficient between each pair of SNP sites in the secondary feature gene pool, retain SNP sites with the Pearson correlation coefficient less than a second set value to obtain a tertiary feature gene pool.

[0020] In a third aspect, the present application further provides an electronic device, characterized in that comprising a memory and a processor, wherein, the memory is configured to store a program; the processor is coupled to the memory and is configured to execute the program stored in the memory to implement steps in the genomic hierarchical feature screening method according to any one of claims 1 to 7 and / or implement steps in the gene screening effect evaluation method according to claim 8.

[0021] The beneficial effects of the present application are: the genomic hierarchical feature screening method provided by the present application firstly acquires original genotype data, then groups the genotype information at each SNP site by chromosomes, fills in the missing genotypes for each chromosome after grouping, and combines all the chromosomes after filling in to obtain a first-level feature gene pool, groups the SNP site genotype information by chromosomes, fills in the missing genotypes after grouping, and combines, perfects the data structure, reduces the influence of data missing, and provides complete and ordered gene data for subsequent screening, calculates the deviation between the combined allele frequency and the expected frequency of two SNP sites on the whole chromosome in the first-level feature gene pool, calculates the linkage disequilibrium coefficient based on the deviation, retains the SNP sites with a linkage disequilibrium coefficient less than a first set value to obtain a second-level feature gene pool, removes redundant sites with high linkage by linkage disequilibrium screening, reduces the data dimension, focuses on SNP sites with independent genetic information, further, calculates the Pearson correlation coefficient between each pair of SNP sites in the second-level feature gene pool, retains the SNP sites with a Pearson correlation coefficient less than a second set value to obtain a third-level feature gene pool, removes the recessive genetic association irrelevant to the physical position through cross-chromosome correlation screening, and significantly improves the genetic independence of the feature pool and the prediction accuracy of the subsequent model. BRIEF DESCRIPTION OF DRAWINGS

[0022] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0023] Figure 1 An embodiment flowchart of a genomic hierarchical feature screening method provided by the present application; Figure 2 An embodiment flowchart of preprocessing of original genotype data provided by the present application; Figure 3 An embodiment flowchart of S102 in the present application; Figure 1 Figure 4 An embodiment flowchart of a gene screening effect evaluation method provided by the present application; Figure 5 An embodiment diagram of the genomic hierarchical feature screening method provided by the present application compared with other feature screening methods in the prior art; Figure 6 An embodiment structure diagram of a genomic hierarchical feature screening device provided by the present application; Figure 7 ​An embodiment structure schematic diagram of the electronic device provided by the present application is shown. DETAILED DESCRIPTION

[0024] The technical solutions in the embodiments of the present application will be clearly and completely described in connection with the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by a person skilled in the art without creative work fall within the protection scope of the present application.

[0025] In the description of the embodiments of the present application, the meaning of "a plurality of" is two or more than two, unless otherwise specified. The association relationship of the associated objects is described by "and / or", which means that there can be three relationships, for example: A and / or B, which can represent the three cases of A existing alone, A and B existing together, and B existing alone.

[0026] The "first", "second", and the like described in the embodiments of the present application are only for the purpose of description, and cannot be understood as indicating or implying the relative importance of the indicated technical features or implicitly indicating the number of the indicated technical features. Therefore, the technical features limited by "first" and "second" can explicitly or implicitly include at least one of the features.

[0027] In this document, the reference to "an embodiment" means that a particular feature, structure, or characteristic described in connection with the embodiment can be included in at least one embodiment of the present application. The appearance of this phrase in various places in the specification does not necessarily all refer to the same embodiment, nor is it necessarily mutually exclusive of other embodiments. It is explicitly and implicitly understood that the embodiments described herein can be combined.

[0028] Before the embodiments are displayed, the following terms are explained.

[0029] Gene (genetic factor): It is the complete nucleotide sequence required to produce a polypeptide chain or functional RNA.

[0030] Chromosome: It refers to a complex substance that can be strongly stained by basic dyes during cell division period, mainly composed of DNA, basic proteins and acidic proteins, and a small amount of RNA.

[0031] SNP site: Single nucleotide polymorphism (SNP) mainly refers to the DNA sequence diversity caused by the variation of a single nucleotide at the genome level, which is numerous and polymorphic.

[0032] Genetic screening: Genetic screening is a technique that detects the genetic information of unborn or newborn babies through genetic technology to avoid major genetic diseases.

[0033] Allele: refers to the gene located at the same position of a pair of homologous chromosomes that controls different forms of the same trait.

[0034] Allele frequency: refers to the proportion of a specific allele at a certain locus in a population, which is used to quantify the genetic diversity or composition of the gene pool of the population.

[0035] Genotype: refers to the total combination of all genes of an individual organism.

[0036] Linkage disequilibrium (LD): refers to the genetic phenomenon of non-random combination of alleles at different loci in a population.

[0037] The present application provides a genomic hierarchical feature screening method, device and electronic equipment, which are described below.

[0038] Figure 1 An embodiment process schematic diagram of a genomic hierarchical feature screening method provided by the present application is shown in Figure 1 The genomic hierarchical feature screening method includes: S101, obtaining original genotype data, the original genotype data including genotype information and sample information at each SNP site of the genome, the genotype information at each SNP site including chromosome number and SNP site number.

[0039] It should be noted that the original genotype data is provided in the binary genotype format of Plink, which is a tool for whole genome association analysis, supporting multiple file formats to store genotype and phenotype data. Common Plink file formats include.ped,.map,.bed,.fam and.bim files.

[0040] It should be noted that the genotype information at each SNP site is stored in a gene matrix, usually in the.bed format, including chromosome number and SNP site number; the sample information includes specific SNP site information and individual information. The specific SNP site information is contained in another text file.bim, with each line containing information of one SNP, including 6 columns: chromosome number, SNP ID (variant name, or. indicates missing), genetic distance (often 0, placeholder), base position (bp unit, physical position), allele 1 (usually reference allele), and allele 2 (usually variant allele). Individual information is stored in.fam files, with each row representing a sample and containing 6 columns: Family ID, Within-family ID, Paternal ID (0 represents unknown), Maternal ID (0 represents unknown), Gender (1 = male, 2 = female, 0 = unknown), and Phenotype value (numeric, -9 or 0 or NA represents missing).

[0041] In some embodiments of the present application, the original genotype data uses the original genotype and corresponding phenotype data (hundred kilograms of age, back fat thickness and teat number) of 2548 Duroc pigs.

[0042] In some embodiments of the present application, the original genotype data uses the original genotype and corresponding phenotype data (chest width and keel length) of 542 Beijing ducks.

[0043] S102, group the genotype information on each SNP site by chromosome, fill in the missing genotype for each chromosome after grouping, and merge all chromosomes after filling to obtain a primary characteristic gene pool.

[0044] It should be noted that: calculating with the whole chromosome can break through the limitation of traditional segmented window, and eliminate the risk of missing screening of high linkage sites across windows.

[0045] In some embodiments of the present application, S103, calculate the deviation between the combined allele frequency and the expected frequency of two SNP sites on the whole chromosome in the primary characteristic gene pool, calculate the linkage disequilibrium coefficient based on the deviation, retain the SNP sites with a linkage disequilibrium coefficient less than a first set value, and obtain a secondary characteristic gene pool.

[0046] In some embodiments of the present application, the first set value is 0.8.

[0047] S104, calculate the Pearson correlation coefficient between each pair of SNP sites in the secondary characteristic gene pool, retain the SNP sites with a Pearson correlation coefficient less than a second set value, and obtain a tertiary characteristic gene pool.

[0048] In some embodiments of the present application, the second set value is 0.9.

[0049] In conclusion, the genome hierarchical feature screening method provided by the application firstly acquires original genotype data, then groups the genotype information at each SNP site by chromosomes, fills in the missing genotypes for each chromosome after grouping, and combines all the chromosomes after filling to obtain a first-level feature gene pool, calculates the deviation between the combined allele frequency and the expected frequency between two SNP sites on the whole chromosome in the first-level feature gene pool, calculates the linkage disequilibrium coefficient based on the deviation, retains the SNP sites with the linkage disequilibrium coefficient less than a first set value to obtain a second-level feature gene pool, since the whole chromosome is calculated, the traditional segmentation window limit can be broken through, and the risk of missing screening of high linkage sites across windows can be eliminated, further, the Pearson correlation coefficient between each pair of SNP sites in the second-level feature gene pool is calculated, and the SNP sites with the Pearson correlation coefficient less than a second set value are retained to obtain a third-level feature gene pool, through cross-chromosome correlation screening, the recessive genetic association irrelevant to the physical position is eliminated, and the genetic independence of the feature pool and the subsequent model prediction accuracy are significantly improved.

[0050] In some embodiments of the application, as shown in Figure 2 Before the step of grouping the genotype information at each SNP site by chromosomes, the step further includes: S201, calculate the minimum allele frequency of each SNP site, and screen out the SNP sites with the minimum allele frequency less than a third set value.

[0051] In some embodiments of the application, the third set value is 0.05.

[0052] S202, perform balance test on each SNP site, and screen out the SNP sites with the significance probability value less than a fourth set value in the balance test.

[0053] In some embodiments of the application, the fourth set value is 0.01.

[0054] S203, calculate the missing rate of each SNP site, and screen out the SNP sites with the missing rate greater than a fifth set value.

[0055] In some embodiments of the application, the fifth set value is 0.01.

[0056] S204, retain all the SNP sites that are not screened out to obtain screened genotype data.

[0057] It should be noted that the screening step before grouping the genotype information at each SNP site by chromosomes is a preprocessing step, which can reduce the influence of sequencing errors and statistical fluctuations.

[0058] In some embodiments of the present application, the original genotype data (hundred kilogram day age, back fat thickness and teat number) of Duroc pigs is used, after a pretreatment step, to obtain a processed Plink binary file, leaving about 11 million SNP sites.

[0059] In some embodiments of the present application, the original genotype data (chest width and keel length) of Beijing ducks is used, after a pretreatment step, to obtain a processed Plink binary file, leaving about 40,000 SNP sites.

[0060] In some embodiments of the present application, the minimum allele frequency of each SNP site in step S201 is calculated by the following formula: (1) (2) In the formula, is the minimum allele frequency of the SNP site, , A, B, etc. represent the frequencies of A, B alleles, respectively, satisfying , N represents the total number of individuals, AA, BB and AB represent the three genotypes that can occur at a certain SNP site, and are the number of AA homozygous genotype, the number of BB homozygous genotype and the number of AB heterozygous genotype, respectively.

[0061] In some embodiments of the present application, the balance test for each SNP site in step S202 includes: calculating the expected frequency of the three genotypes at each SNP site, the expected frequency of the AA homozygous genotype is , the expected frequency of the BB homozygous genotype is , and the expected frequency of the AB heterozygous genotype is ; calculating the expected number of the three genotypes at each SNP site, the expected number is the expected frequency multiplied by the total number of individuals; calculating the statistical quantity of the chi-square test of the actual observed number and the expected number of the three genotypes χ 2 , by the following formula: (3) In the formula, , and are the expected numbers of the three genotypes; calculating the degrees of freedom of the chi-square test, the degrees of freedom is the difference between the number of genotypes and the number of alleles; The significance probability value of the chi-square test is obtained by referring to the chi-square distribution critical value table through the chi-square test statistic and degrees of freedom.

[0062] In some embodiments of the present application, in step S203, the deletion rate of each SNP site is calculated as the quotient of the number of deletion samples of the SNP site and the total number of samples.

[0063] In some embodiments of the present application, as shown in step S102, the sample information includes a phenotype value; the genotype information on each SNP site is grouped by chromosome, the missing genotype of each grouped chromosome is filled, and all the filled chromosomes are merged by chromosome to obtain a primary feature gene pool, which includes: Figure 3 S301, grouping each SNP site in the screened genotype data according to the chromosome number and saving it as a separate file for each chromosome; S302, estimating the possible genotype of the missing SNP site in each chromosome through the linkage block in the screened genotype data, and filling the genotype of the missing SNP site in the chromosome; S303, merging the filled chromosomes to obtain a complete genotype matrix, and performing whole genome association analysis according to the phenotype value, screening the SNP sites with a significance probability value less than a sixth set value in the whole genome association analysis, and obtaining a primary feature gene pool; It should be noted that the linkage block refers to the SNP sites closely arranged in the chromosome.

[0064] It should be noted that grouping each SNP site in the screened genotype data according to the chromosome number can use the Plink tool to convert the binary file to a vcf file, and then use the grep tool to group the vcf primary feature pool converted by Plink according to the chromosome number and save it as a separate file for each chromosome, which is convenient for parallel processing.

[0065] It should be noted that the genotype of the missing SNP site in the chromosome can be filled by using a hidden Markov model (HMM) to infer the linkage block in the population, thereby estimating the possible genotype state of the missing site, and then using the Beagle software to fill the missing genotype of each chromosome in the partitioned primary feature pool in parallel with 4 threads, i.e. setting nthreads=4.

[0066] In some embodiments of the present application, the original genotype data (hundred-kilogram day age, back fat thickness and teat number) of Duroc pigs is used to fill the genotype of the missing SNP site in the chromosome, and 18 vcf format chromosome files are obtained.

[0067] ​In some embodiments of the present application, the raw genotype data (chest width and keel length) of Beijing duck is used to fill in the genotype of the missing SNP sites in the chromosome, and 29 chromosome files in vcf format are obtained.

[0068] In some embodiments of the present application, the sixth set value is 0.001.

[0069] In some embodiments of the present application, the raw genotype data (hundred kilogram day age, back fat thickness and teat number) of Duroc pig is used to screen SNP sites with a significance probability value less than the sixth set value in the genome-wide association analysis, and a primary characteristic gene pool is obtained. After processing, the primary characteristic pools of hundred kilogram day age, back fat thickness and teat number are left with 332,870, 412,631 and 478,646 SNP sites, respectively.

[0070] In some embodiments of the present application, the raw genotype data (chest width and keel length) of Beijing duck is used to screen SNP sites with a significance probability value less than the sixth set value in the genome-wide association analysis, and a primary characteristic gene pool is obtained. After processing, the primary characteristic pools of chest width and keel length are left with 2475 and 2255 SNP sites, respectively.

[0071] In some embodiments of the present application, the deviation between the joint allele frequency and the expected frequency of two SNP sites on the whole chromosome in the primary characteristic gene pool is calculated, and the linkage disequilibrium coefficient is calculated based on the deviation, including: The deviation between the joint allele frequency and the expected frequency is calculated , which is calculated by the following formula: (4) In the formula, is the frequency of the simultaneous occurrence of haploid genotypes A and B, is the frequency of allele A, is the frequency of allele B; The deviation is used to calculate the linkage disequilibrium coefficient between two SNP sites on each chromosome , which is calculated by the following formula: (5) In the formula, is the frequency of non-A allele of the first SNP site, is the frequency of non-B allele of the second SNP site.

[0072] In some embodiments of the present invention, this step can use the --r2 command of Plink v1.90b6.21 software and set the --ld-window 999999 parameter to calculate the whole chromosome, and screen out sites with low linkage disequilibrium coefficients from the primary feature pool to form a secondary feature pool, wherein the chromosome with the higher minimum allele frequency of the SNP site is retained.

[0073] It should be noted that: calculating on the whole chromosome can ensure that no long-distance linkage disequilibrium is missed within the chromosome, and can also avoid the difference in results caused by the window parameter (--ld-window), minimize the number of linkage disequilibrium sites on each chromosome, and remove redundant genes on each chromosome.

[0074] In some embodiments of the present invention, after performing linkage disequilibrium screening using the original genotype data (age per 100 kg, backfat thickness, and number of teats) of Duroc pigs, the secondary characteristic gene pools for age per 100 kg, backfat thickness, and number of teats have 7158, 8246, and 12219 SNP loci remaining, respectively.

[0075] In some embodiments of the present invention, after using the original genotype data (breast width and keel length) of Beijing ducks, the secondary characteristic gene pools for breast width and keel length obtained by linkage disequilibrium screening have 2209 and 2068 SNP loci remaining, respectively.

[0076] In some embodiments of the present invention, calculating the Pearson correlation coefficient between each pair of SNP sites in the secondary feature gene pool includes: The Pearson correlation coefficient between each pair of SNP loci is calculated using the following formula: (6) In the formula, and For the first i The gene coding value of each individual at the corresponding SNP locus. and These are the average gene coding values ​​of two SNP loci. N This represents the total number of individuals.

[0077] It should be noted that calculating the Pearson correlation coefficient between each pair of SNP loci, and calculating the correlation between loci, and filtering according to the set values, can minimize physical linkage disequilibrium.

[0078] In some embodiments of the present invention, using the original genotype data of Duroc pigs (age per 100 kg, backfat thickness, and number of teats), after calculating the Pearson correlation coefficient between each pair of SNP loci and screening, a three-level characteristic gene pool is obtained. The three-level characteristic gene pools for age per 100 kg, backfat thickness, and number of teats have 1111, 1056, and 1393 SNP loci remaining, respectively.

[0079] In some embodiments of the present invention, using the original genotype data (breast width and keel length) of Beijing ducks, after calculating the Pearson correlation coefficient between each pair of SNP loci and screening, a three-level characteristic gene pool is obtained. The three-level characteristic gene pools for breast width and keel length have 1933 and 1812 SNP loci remaining, respectively.

[0080] To evaluate the effectiveness of genome-level feature screening methods, such as Figure 4 As shown, the present invention also provides a method for evaluating the effect of gene screening, comprising: S401. A three-level feature gene pool is determined based on a genome-level feature screening method.

[0081] S402. Based on the mapping relationship between the three-level feature gene pool and the gene pool-phenotype value, predict the unknown phenotypic values ​​of all SNP sites to obtain all predicted phenotypic values.

[0082] It should be noted that the gene pool-phenotypic value mapping relationship is used to characterize the correspondence between the tertiary feature gene pool and the phenotypic values ​​corresponding to SNP sites in the tertiary feature gene pool.

[0083] S403. Calculate the correlation coefficient between the predicted total phenotypic values ​​and the actual total phenotypic values, and evaluate the gene screening effect using the correlation coefficient.

[0084] In some embodiments of the present invention, the correlation coefficient between the predicted total phenotypic values ​​and the actual total phenotypic values ​​can be calculated using the Pearson correlation coefficient, the formula for which the Pearson correlation coefficient is calculated is as follows: (7) In the formula, For the first i Predicted phenotypic values ​​for each individual For the first i The true phenotypic value of each individual To predict the phenotypic mean, This represents the true phenotypic mean.

[0085] It should be noted that the breeding value is estimated using existing statistical models to predict the target trait of individuals with unknown phenotypes and genotypes. The input is a matrix of known phenotypes and genotypes, and the output is the phenotypes of all individuals. The unknown is deduced from the known, and a complete judgment of the breeding genes is made.

[0086] In some embodiments of the present application, the existing statistical models are GBLUP, RRGBLUP, BayesLasso, BayesA, ByaesB and BayesC.

[0087] It should be noted that: GBLUP is a genomic best linear unbiased prediction model, RR-GBLUP is a ridge regression genomic best linear unbiased prediction model, BayesLasso is a Bayesian lasso regression model, BayesA is a Bayesian A model, ByaesB is a Bayesian B model, and BayesC is a Bayesian C model.

[0088] It should be noted that: in order to show the effectiveness of the genomic hierarchical feature screening method, when using the existing statistical model to estimate the breeding value, other screening methods will also be used to estimate the breeding value by using the existing statistical model, so as to compare the screening methods.

[0089] In some embodiments of the present application, as shown in Figure 5 The examples of different screening methods are A-F in the figure, and A-F are respectively: A: LD screening, only single LD linkage disequilibrium screening is performed; B: LD+GWAS, LD screening retains important SNP sites analyzed by genome-wide association analysis (GWAS); C: LD+Epi, LD screening retains high-effect SNP sites calculated by SNP sites in Epi; D: LD+PCA, further principal component analysis (PCA) is performed after LD screening; E: LD+outPCA, the excluded sites are added to the feature pool after principal component analysis; F: LD+PCC, the genomic hierarchical feature screening method of the present application.

[0090] Among them, the A method only analyzes linkage disequilibrium within the chromosome, eliminates high linkage (redundancy) sites, can quickly reduce the number of redundant sites from homologous chromosome regions (i.e. within the chromosome), reduce the data dimension, but cannot identify the association between sites across chromosomes (implicit redundancy); the B method can preferentially retain sites with significant additive effects on traits, ensuring that core functional sites are not lost, and may improve the prediction ability of major genes, but the effect is highly dependent on the accuracy and effectiveness of the previous GWAS analysis, GWAS itself is difficult to detect microgenes, epistatic effects, etc. Non-additive effect sites may be eliminated; the C method relies on significant sites, can capture interaction effects beyond additive effects, dominant effects, etc. More consistent with the genetic architecture of complex traits, but the sample size required for epistatic effect detection is extremely large, the statistical power is very low in reality, the results are unstable, and false positive interaction signals are extremely easy to produce, resulting in unreliable retained sites; the D method principal component analysis can extract the maximum variance direction, and a few principal components represent most of the genetic variation, greatly improving the calculation efficiency, but it may filter out rare variations that have a weak but critical impact on specific traits; the E method adds the feature pool after principal component analysis of the screened sites, which can prevent some important sites from being removed, but the removed sites are added back in a transformed form, which may reintroduce redundancy or noise, which violates the original intention of dimension reduction; F is the method of the present application, which can reduce the dimension of the chromosome, and through local and global dimension reduction, double redundancy reduction, the dimension reduction is more thorough, although two rounds of calculation are required, but thanks to the previous dimension reduction, the total is still within an acceptable range.

[0091] It should be noted that "dimension reduction" means reducing the number of SNP sites used to construct a prediction model while retaining as much effective information as possible, resulting in a high-quality, low-redundancy feature set; "microgene" refers to multiple genes that have a small phenotypic effect on quantitative traits; "epistatic effect" refers to the effect caused by the interaction of genes at different sites; "additive effect" refers to the cumulative effect of multiple alleles on a trait, with each allele contributing to the trait; "dominant effect" describes the deviation of the phenotype of a heterozygote (such as Aa) from the two homozygotes.

[0092] In some embodiments of the present application, the original genotype data uses the original genotypes (hundred kilogram day age, back fat thickness and teat number) of 2548 Duroc pigs for gene screening, and the number of SNP sites remaining after different screening methods is shown in Table 1.

[0093] Table 1: Number of SNP sites remaining after different gene screening methods (Duroc pigs)

[0094] In some embodiments of the present application, the original genotype data uses the original genotype (breast width and keel length) of 542 Beijing ducks for gene screening, and the number of SNP sites remaining after screening by different screening methods is shown in Table 2.

[0095] Table 2: Number of SNP sites remaining after screening by different gene screening methods (Beijing ducks)

[0096] It should be noted that: in the screening process, the statistics of the number of remaining SNP sites is only the record of experimental data, and the number cannot directly reflect the good or bad of the screening effect.

[0097] In some embodiments of the present application, the original genotype data uses the original genotype and corresponding phenotype data (hundred kilograms of age, back fat thickness and teat number) of 2548 Duroc pigs, and the corresponding three-level characteristic gene pool is obtained by the genomic hierarchical characteristic screening method. According to the mapping relationship between the three-level characteristic gene pool and the gene pool-phenotype value, the unknown phenotype values in all SNP sites are predicted to obtain the predicted all phenotype values, and the Pearson correlation coefficient between the predicted all phenotype values and the actual all phenotype values is calculated. In some embodiments of the present application, the original genotype data uses the original genotype and corresponding phenotype data (breast width and keel length) of 542 Beijing ducks, and the corresponding three-level characteristic gene pool is obtained by the genomic hierarchical characteristic screening method. According to the mapping relationship between the three-level characteristic gene pool and the gene pool-phenotype value, the unknown phenotype values in all SNP sites are predicted to obtain the predicted all phenotype values, and the Pearson correlation coefficient between the predicted all phenotype values and the actual all phenotype values is calculated.

[0098] In order to better reflect the effect of the genomic hierarchical characteristic screening method provided by the present application and other characteristic screening methods in the prior art after prediction, the Pearson correlation coefficients of the predicted phenotype values and the actual phenotype values of all statistical prediction methods are averaged.

[0099] In some embodiments of the present application, the Pearson correlation coefficient results calculated by the screening method and the corresponding prediction method of the original genotype data of 2548 Duroc pigs and the corresponding phenotype data (hundred kilograms of age, back fat thickness and teat number) are shown in Table 3.

[0100] Table 3: Average Pearson correlation coefficient of prediction results of Duroc pig original gene data after gene screening

[0101] It should be noted that: STD in the table is the standard deviation, which is used to describe the aggregation degree of data points around the mean value.

[0102] It should be noted that: it should be noted that the screening methods in the table are respectively: Orgin: represents the first-level feature pool after screening; LD: only LD linkage disequilibrium screening alone; LD_GWASp10-10: LD screening retains important SNP sites for genome-wide association analysis (GWAS), and the data after GWAS represents the threshold of P value, p10-10 represents P<1×10 -10 ; LD_GWASp10-12: LD screening retains important SNP sites for genome-wide association analysis (GWAS), and the data after GWAS represents the threshold of P value, p10-12 represents P<1×10 -12 ; LD_PCA: further PCA after LD screening; LD_PCC: the genome-level feature screening method of the application; LD_Epi: LD screening retains high-effect SNP sites calculated from SNP sites in EPI.

[0103] Among them, the gene screening method (LD_PCC) of the application on the original genotype data of the Duroc pig's hundred kilogram day age, back fat thickness and nipple number, the average of the Pearson correlation coefficient of the prediction result after gene screening is the highest among all the compared gene screening methods. The higher the Pearson correlation coefficient here means that the individual with high predicted value also has high true value, and the individual with low predicted value also has low true value, which can indicate that the average accuracy of the prediction of the gene screening method of the application is higher, that is, the gene screening method of the application is more effective.

[0104] In some embodiments of the application, the original genotype data of 542 Beijing ducks is used to calculate the Pearson correlation coefficient of the screening method and the corresponding prediction method of the original genotype and corresponding phenotype data (breast width and keel length). The results are shown in Table 4.

[0105] Table 4: The average of the Pearson correlation coefficient of the prediction result after the original gene data of Beijing duck is screened by gene

[0106] It should be noted that: it should be noted that the screening methods in the table are respectively: Orgin: represents the first-level feature pool after screening; LD: only LD linkage disequilibrium screening alone; LD_GWASp1: LD screening retains important SNP sites for genome-wide association analysis (GWAS), and the data after GWAS represents the threshold of P value, p1 represents P < 0.01; LD_GWASp2: LD screening retains important SNP sites for genome-wide association analysis (GWAS), and the data after GWAS represents the threshold of P value, p2 represents P < 0.02; LD_PCA: further PCA after LD screening; LD_PCC: genome hierarchical feature screening method of the application; LD_Epi: high-effect SNP sites calculated by SNP sites in EPI during LD screening.

[0107] Among them, the gene screening method (LD_PCC) of the application obtains the highest value (keel length) and the second highest value (breast width) in the average Pearson correlation coefficient of the prediction results after gene screening on the original genotype data of the breast width and keel length of Beijing ducks among all the compared gene screening methods. The higher the Pearson correlation coefficient here means that the individual with high predicted value also has high true value, and the individual with low predicted value also has low true value, which can indicate that the average prediction accuracy of the gene screening method of the application is higher, that is, the gene screening method of the application is more effective.

[0108] In order to better implement the genome hierarchical feature screening method in the embodiment of the application, on the basis of the genome hierarchical feature screening method, as shown in Figure 6 The embodiment of the application further provides a genome hierarchical feature screening device, and the genome hierarchical feature screening device 600 comprises: An original genotype data acquisition unit 601 is configured to acquire original genotype data, wherein the original genotype data comprises genotype information and sample information of each SNP site in a genome, and the genotype information of each SNP site comprises a chromosome number and a SNP site number. A primary feature gene screening unit 602 is configured to group the genotype information of each SNP site by chromosomes, fill in the missing genotypes of each chromosome after grouping, and combine all the chromosomes after filling to obtain a primary feature gene pool. A secondary feature gene screening unit 603 is configured to calculate the deviation between the joint allele frequency and the expected frequency of two SNP sites on an entire chromosome in the primary feature gene pool, calculate the linkage disequilibrium coefficient based on the deviation, retain the SNP sites with the linkage disequilibrium coefficient less than a first set value, and obtain a secondary feature gene pool. The tertiary feature gene screening unit 604 is used to calculate the Pearson correlation coefficient between each pair of SNP sites in the secondary feature gene pool, and retain SNP sites with a Pearson correlation coefficient less than a second set value to obtain the tertiary feature gene pool.

[0109] The genome hierarchical feature screening device 600 provided in the above embodiments can realize the technical solutions described in the above genome hierarchical feature screening method embodiments. The specific implementation principles of each module or unit can be found in the corresponding content in the above genome hierarchical feature screening method embodiments, and will not be repeated here.

[0110] like Figure 7 As shown, the present invention also provides an electronic device 700. The electronic device 700 includes a processor 701, a memory 702, and a display 703. Figure 7 Only some components of the electronic device 700 are shown, but it should be understood that it is not required to implement all of the components shown, and more or fewer components may be implemented instead.

[0111] In some embodiments, processor 701 may be a central processing unit (CPU), microprocessor, or other data processing chip, used to run program code stored in memory 702 or process data, such as the genome hierarchical feature screening method and gene screening effect evaluation method in this invention.

[0112] In some embodiments, processor 701 may be a single server or a group of servers. The server group may be centralized or distributed. In some embodiments, processor 701 may be local or remote. In some embodiments, processor 701 may be implemented on a cloud platform. In one embodiment, the cloud platform may include a private cloud, public cloud, hybrid cloud, community cloud, distributed cloud, intranet, multi-cloud, etc., or any combination thereof.

[0113] In some embodiments, memory 702 may be an internal storage unit of electronic device 700, such as a hard disk or memory of electronic device 700. In other embodiments, memory 702 may also be an external storage device of electronic device 700, such as a plug-in hard disk, smart media card (SMC), secure digital (SD) card, flash card, etc. equipped on electronic device 700.

[0114] Furthermore, the memory 702 may include both internal storage units of the electronic device 700 and external storage devices. The memory 702 is used to store application software and various types of data installed on the electronic device 700.

[0115] The display 703 can be, in some embodiments, an LED display, a liquid crystal display, a touch liquid crystal display, an OLED (Organic Light-Emitting Diode) touch, etc. The display 703 is used to display information of the electronic device 700 and to display a visualized user interface. The components 701-703 of the electronic device 700 communicate with each other through a system bus.

[0116] In an embodiment, when the processor 701 executes the genomic hierarchical feature screening program in the memory 702, the following steps can be implemented: Obtaining raw genotype data, the raw genotype data including genotype information at each SNP site on a genome and sample information, the genotype information at each SNP site including a chromosome number and a SNP site number; Grouping the genotype information at each SNP site by chromosomes, filling in missing genotypes for each chromosome after grouping, and merging all chromosomes after filling in to obtain a primary feature gene pool; Calculating a deviation between a joint allele frequency and an expected frequency between two SNP sites on an entire chromosome in the primary feature gene pool, calculating a linkage disequilibrium coefficient based on the deviation, and retaining SNP sites with linkage disequilibrium coefficients less than a first set value to obtain a secondary feature gene pool; Calculating a Pearson correlation coefficient between each pair of SNP sites in the secondary feature gene pool, retaining SNP sites with Pearson correlation coefficients less than a second set value to obtain a tertiary feature gene pool.

[0117] In an embodiment, when the processor 701 executes the gene screening effect evaluation program in the memory 702, the following steps can be implemented: Determining a tertiary feature gene pool based on the genomic hierarchical feature screening method; Predicting unknown phenotype values in all SNP sites according to a mapping relationship between the tertiary feature gene pool and a gene pool-phenotype value, to obtain predicted all phenotype values; the gene pool-phenotype value mapping relationship is used to represent a corresponding relationship between the tertiary feature gene pool and phenotype values corresponding to SNP sites in the tertiary feature gene pool; Calculating a correlation coefficient between the predicted all phenotype values and actual all phenotype values, and evaluating a gene screening effect through the correlation coefficient.

[0118] It should be understood that, in addition to the above functions, the processor 701 can also implement other functions when executing the genomic hierarchical feature screening program and the gene screening effect evaluation program in the memory 702. For details, refer to the descriptions of the corresponding method embodiments above.

[0119] Further, the embodiments of the present application do not make specific limitation on the type of the electronic device 700 mentioned above, and the electronic device 700 can be a portable electronic device such as a mobile phone, a tablet computer, a personal digital assistant (PDA), a wearable device, a laptop, etc. Exemplary embodiments of the portable electronic device include, but are not limited to, a portable electronic device running an IOS, an android, a microsoft or other operating system. The portable electronic device mentioned above can also be other portable electronic devices such as a laptop having a touch-sensitive surface (e.g., a touch panel), etc. It should also be understood that in some other embodiments of the present application, the electronic device 700 can not be a portable electronic device, but a desktop computer having a touch-sensitive surface (e.g., a touch panel).

Claims

1. A genome-wide feature screening method, characterized by, The method comprises the following steps: obtaining original genotype data, wherein the original genotype data comprises genotype information of each SNP site on a genome and sample information, and the genotype information of each SNP site comprises a chromosome number and a SNP site number; grouping the genotype information of each SNP site by chromosome, filling missing genotypes for each chromosome after grouping, and merging all chromosomes after filling to obtain a first characteristic gene pool; calculating a deviation between a combined allele frequency and an expected frequency between two SNP sites on an entire chromosome in the first characteristic gene pool, calculating a linkage disequilibrium coefficient based on the deviation, and retaining SNP sites with a linkage disequilibrium coefficient less than a first set value to obtain a second characteristic gene pool; calculating a Pearson correlation coefficient between each pair of SNP sites in the second characteristic gene pool, and retaining SNP sites with a Pearson correlation coefficient less than a second set value to obtain a third characteristic gene pool.

2. The genomic hierarchical feature screening method according to claim 1, characterized by, Before grouping the genotype information of each SNP site by chromosome, the method further comprises the following steps: calculating a minimum allele frequency of each SNP site, and screening out SNP sites with a minimum allele frequency less than a third set value; performing a balance test on each SNP site, and screening out SNP sites with a significant probability value less than a fourth set value in the balance test; calculating a missing rate of each SNP site, and screening out SNP sites with a missing rate greater than a fifth set value; retaining all SNP sites that are not screened out to obtain screened genotype data.

3. The genomic hierarchical feature screening method according to claim 2, wherein, The minimum allele frequency of each SNP site is obtained by the following calculation formula: wherein, is the minimum allele frequency at the SNP site, , respectively represent the frequencies of A, B alleles, satisfying , N represents the total number of individuals, AA, BB and AB represent the three possible genotypes at a SNP site, and are respectively the number of AA homozygous genotype, the number of BB homozygous genotype and the number of AB heterozygous genotype.

4. The genomic hierarchical feature screening method according to claim 2, wherein, The balance test on each SNP site comprises the following steps: the expected frequency of the AA homozygous genotype is the expected frequency of the BB homozygous genotype is the expected frequency of the AB heterozygous genotype is ; calculating an expected number of the three genotypes on each SNP site, wherein the expected number is an expected frequency multiplied by a total number of individuals; calculating a statistic for a chi-square test of the actual observed number of three genotypes and the expected number calculating a degree of freedom of the chi-square test, wherein the degree of freedom is a difference between a genotype number and an allele number; 2 is obtained by the following calculation formula: wherein , and is the expected number of the three genotypes; obtaining a significant probability value of the chi-square test by querying a chi-square distribution critical value table according to a statistic value of the chi-square test and the degree of freedom. The sample information comprises a phenotype value; the grouping of the genotype information of each SNP site by chromosome, the filling of missing genotypes for each chromosome after grouping, and the merging of all chromosomes after filling to obtain a first characteristic gene pool comprise the following steps:

5. The genomic hierarchical feature screening method according to claim 2, wherein, grouping each SNP site in the screened genotype data according to a chromosome number, and saving each chromosome as a separate file; estimating possible genotypes of missing SNP sites in each chromosome through a linkage block in the screened genotype data, and filling genotypes of missing SNP sites in the chromosome; merging the filled chromosomes to obtain a complete genotype matrix, performing whole genome association analysis according to the phenotype value, screening SNP sites with a significant probability value less than a sixth set value in the whole genome association analysis, and obtaining the first characteristic gene pool; the linkage block refers to closely arranged SNP sites in a chromosome. ​ 6. The genomic hierarchical feature screening method according to claim 1, wherein, the computing the deviation between the joint allele frequency and the expected frequency between two SNP loci on the whole chromosome in the primary feature gene pool, and the computing the linkage disequilibrium coefficient based on the deviation, comprise: The deviation between the combined allele frequency and the expected frequency is calculated is calculated by the following formula: wherein is the frequency of the simultaneous occurrence of haploid genotypes A and B, is the frequency of allele A, is the frequency of allele B; by said deviation , the linkage disequilibrium coefficient between two SNP sites on each chromosome is calculated , by the following formula: wherein is the frequency of non-A alleles at the first SNP site, is the frequency of non-B alleles at the second SNP site.

7. The genomic hierarchical feature screening method according to claim 1, wherein, the computing the Pearson correlation coefficient between each pair of SNP loci in the secondary feature gene pool, comprises: the Pearson correlation coefficient between each pair of SNP loci is calculated by the following formula: wherein and is the gene code value of the i-th individual at the corresponding SNP site, i is the gene code value of the i-th individual at the corresponding SNP site, and are the average gene code values of the two SNP sites, respectively, N is the total number of individuals.

8. A method for evaluating the effect of a gene screening, characterized by, comprise: determining a tertiary feature gene pool based on the genomic hierarchical feature screening method according to any one of claims 1-7; predicting unknown phenotype values in all SNP loci according to the mapping relationship between the tertiary feature gene pool and the gene pool-phenotype value, to obtain predicted all phenotype values; the gene pool-phenotype value mapping relationship is used to represent the correspondence between the tertiary feature gene pool and the phenotype values corresponding to the SNP loci in the tertiary feature gene pool; the correlation coefficient between the predicted all phenotype values and the actual all phenotype values is calculated, and the gene screening effect is evaluated by the correlation coefficient.

9. A genomic hierarchical feature screening device, characterized by, comprise: an original genotype data acquisition unit for acquiring original genotype data, the original genotype data comprising genotype information and sample information of each SNP locus in the genome, the genotype information of each SNP locus comprising chromosome number and SNP locus number; a primary feature gene screening unit for grouping the genotype information of each SNP locus by chromosome, filling the missing genotype for each chromosome after grouping, and merging all chromosomes after filling to obtain a primary feature gene pool; a secondary feature gene screening unit for calculating the deviation between the joint allele frequency and the expected frequency between two SNP loci on the whole chromosome in the primary feature gene pool, and calculating the linkage disequilibrium coefficient based on the deviation, and retaining SNP loci with linkage disequilibrium coefficient less than a first set value to obtain a secondary feature gene pool; a tertiary feature gene screening unit for calculating the Pearson correlation coefficient between each pair of SNP loci in the secondary feature gene pool, and retaining SNP loci with Pearson correlation coefficient less than a second set value to obtain a tertiary feature gene pool.

10. An electronic device, comprising: comprise a memory and a processor, wherein the memory is used to store programs; the processor is coupled with the memory, and is used to execute the programs stored in the memory to realize the steps in the genomic hierarchical feature screening method according to any one of claims 1-7, and / or realize the steps in the gene screening effect evaluation method according to claim 8.

Images