High genetic information SNP locus set screening method and device, medium and computer program product

Abstract

The invention discloses a high genetic information SNP (Single Nucleotide Polymorphism) site set screening method and equipment, a medium and a computer program product. According to the screening method, similarity genome intervals are divided through an LD threshold value and an optimized K-means model, and representative SNP loci in the divided intervals are screened out by applying a variance analysis technology; after screening is completed, a high genetic information SNP site set is output. The invention provides equipment, a medium and a computer program product capable of realizing the screening method. The screening method is wide in application range, directly selects the SNP site set with the largest genomic genetic information amount, is high in interpolation accuracy, and can effectively reduce the design and production cost while greatly improving the effectiveness of chip sites.

Description

Technical Field

[0001] The present invention relates to a SNP site screening method, and in particular to a high genetic information SNP site set screening method, equipment, medium and computer program product. Background Art

[0002] Single nucleotide polymorphism (SNP) is the most common form of genetic variation in the genome, and is closely related to individual traits, disease susceptibility, population structure analysis, and evolutionary research. Screening SNP sites with high genetic information is crucial for precision medicine, animal and plant genetic breeding, and biodiversity research.

[0003] Traditional SNP site screening methods are based on whole genome data calculations, but are easily interfered by redundant information, have low screening efficiency, and the selected sites may be insufficiently stable or contain low information in practical applications, making it difficult to meet the needs of efficient analysis. Especially in animal genetic breeding, in order to ensure that the chip has high applicability and reliability, SNP sites that are widely distributed among breeds or have known functions are usually selected. Since the differences in genomic haplotype patterns among breeds are not taken into account, these chips often have problems such as limited representativeness of marker sites, failure of some sites, and difficulty in detecting breed-specific genetic information in the application of specific populations. Summary of the invention

[0004] Purpose of the invention: The purpose of the present invention is to provide a method for screening a set of SNP sites with high genetic information; the second purpose is to provide equipment, media and computer program products that can implement the screening method.

[0005] Technical solution: A method for screening a set of high genetic information SNP sites according to the present invention comprises:

[0006] (1) Divide haploid domains based on LD threshold;

[0007] (2) The SNP sites not included in the monomer domain interval were divided by K-means model;

[0008] (3) Using variance analysis technology, representative SNP sites within the similar genomic intervals divided in the previous step were screened out;

[0009] (4) After the screening is completed, a set of SNP sites with high genetic information is output.

[0010] Preferably, in the method for screening a set of high genetic information SNP sites of the present invention, dividing the haploid domain based on the LD threshold in step 1 comprises:

[0011] (11) Calculate the LD value between any two SNP sites in the genome of each individual;

[0012] (12) The calculation results are arranged in order according to the position of the SNP sites in the genome to form a symmetric matrix;

[0013] (13) Setting the LD threshold, taking the first SNP site of the first chromosome as the starting point, merging sites to the right, and calculating the minimum value of all elements in the square matrix formed by each starting SNP and the newly merged SNP in the aforementioned matrix until the minimum value of all elements in the square matrix formed by the starting SNP and the newly merged SNP is less than the threshold;

[0014] (14) Remove the most recently merged SNP site, and use the starting SNP site and the last merged SNP site as the interval. If there are two or more SNP sites in the interval, the interval is regarded as a haploid region;

[0015] (15) Using the SNP closest to the right end of the newly divided haplotype block as a new starting point, repeat steps 13-14 until the genome is completely divided.

[0016] Preferably, in the method for dividing haploid domains based on LD threshold value of the present invention, the LD threshold value in step 13 is a specific value between 0 and 1.

[0017] Preferably, in the method for screening a set of high genetic information SNP sites of the present invention, in step 2, the SNP sites not included in the monomer domain interval are divided by an optimized K-means model, comprising:

[0018] (21) SNP sites that were not included in the monomer domain interval were integrated in order of position;

[0019] (22) The amount of genetic information of the above SNP site is calculated by formula 221, the conditional entropy under the random variable site is calculated by formula 222, and the normalized information gain is calculated by formula 223:

[0020]

[0021] Among them, H(SNP x )、H(SNP x |SNP y ) and SU(SNP x ,SNP y ) are SNPs x Genetic information, SNP x In the random variable SNP y Conditional entropy and normalized information gain under P(SNP x =i) is SNP x The probability of genotype i;

[0022] (23) The Euclidean distance between the above SNP sites is calculated using the modified Euclidean distance metric formula 231, and the density distance value of all SNP sites is calculated using formula 232:

[0023]

[0024] Among them, SNP xt and SNP yt is the genotype value of the corresponding t-th individual; SU(SNP x ,SNP y ) is the normalized information gain; e is the natural logarithm; r 1 is the effective density radius, which is usually replaced by the field radius; T is the total number of data points within the radius;

[0025] (24) According to the elbow rule, determine the optimal number of clusters K 1 , select step 23 to calculate the first K with the largest density value 1 The SNP sites are used as the initial cluster centers. Based on the improved Euclidean distance measurement formula, the distances between all SNP sites not included in the monomer domain interval and the cluster center vectors are calculated again, and they are classified into the category of the nearest cluster center.

[0026] (25) Calculate K 1 The sum of squared errors of the Euclidean distances between SNP sites in each cluster is calculated. The SNP site with the smallest absolute value of the sum of squared errors of all SNP sites in each cluster is used as the new cluster center. The distances between all SNP sites in the genome and the new K 1 The Euclidean distances of cluster centers are re-divided into K 2 Clusters. Repeat this step until the cluster center no longer changes.

[0027] Preferably, in the method for screening a set of high genetic information SNP sites of the present invention, the step 3 of screening sites by using variance analysis technology comprises:

[0028] (31) For each divided interval in step 1-2, the SNP site data is converted into 0, 1 or 2, and the genotype scaling value of each SNP site is calculated by formula 311:

[0029]

[0030] in, and θ(SNP α ) are SNPs α The mean and standard deviation of the genotype values ​​of the locus in the population;

[0031] (32) Based on the genotype scaling value, the covariance matrix between the SNP sites in each divided interval is calculated, and the correction value of the amount of genomic information explained by each SNP site after removing the interference of other SNP sites in the interval is calculated by formula 321:

[0032]

[0033] Among them, V δδ 、V αα and V δα SNP δ and SNP α The variance and covariance between them, N is the number of SNP sites in the interval

[0034] (33) Select the SNP site with the largest correction value of the amount of genome information explained in each genome division interval as the representative site of the interval, and use formula 331 to calculate the (co)variance correction matrix between the remaining SNP sites in the interval. After correction, screen the next round of candidate high genetic information SNP sites in the interval;

[0035]

[0036] Among them, V βγ Select SNPs for this interval ξ The scaling matrix before the site; V βξ 、V γξ and V ξξ Other SNP sites and SNP in the interval are ξ The (co)variance size, N is the number of SNP sites in this interval;

[0037] (34) After each round of screening of candidate high genetic information SNP sites, the total proportion of genomic information in each genome division interval is calculated by formula 341:

[0038]

[0039] Among them, θ 0 2 is the total variance of the matrix before selecting the SNP site, θ t 2 is the residual variance of the matrix after selecting SNP sites, t is the number of selected sites, and for each divided interval, a threshold w is set. When the total proportion of genomic information explained by the SNP sites selected in this interval reaches the threshold, it stops and slides to the next genome division interval for screening.

[0040] Preferably, in the method for selecting sites using the variance analysis technology of the present invention, in step 32, one or more SNP sites with the largest genomic information correction value in each genomic division interval are selected.

[0041] Preferably, in the method for screening sites using the variance analysis technology of the present invention, the threshold value w in step 34 is a specific value between 0.6 and 1.

[0042] The computer device of the present invention comprises a memory, a processor and a computer program stored in the memory, and the processor executes the steps of the high genetic information SNP site set screening method.

[0043] The computer-readable storage medium of the present invention stores a computer program / instruction thereon, which, when executed by a processor, implements the steps of the method for screening a set of high genetic information SNP sites.

[0044] The computer program product of the present invention includes a computer program / instruction, which, when executed by a processor, implements the steps of a method for screening a set of high genetic information SNP sites.

[0045] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: 1. It has a wide range of applications and can be applied as long as there is genome sequencing or SNP site information; 2. It directly selects the SNP site set with the largest amount of genomic genetic information, which greatly improves the effectiveness of breeding chip sites and reduces chip design and production costs; 3. The screened SNP set has good interpolation characteristics, and the interpolation accuracy of 140K sites reaches 0.86±0.02, which is significantly higher than the accuracy of 100K sites commonly used on the market. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 Provide a technical roadmap for the screening method of high genetic information SNP loci sets;

[0047] Figure 2 This is a comparison chart of the interpolation accuracy of the 10K site set;

[0048] Figure 3 This is a comparison chart of the interpolation accuracy of the 140K site set;

[0049] Figure 4 A schematic diagram of the structure of a computer device. DETAILED DESCRIPTION

[0050] The technical solution of the present invention is further described below.

[0051] Example 1: Screening method for high genetic information SNP loci set for screening high genetic information loci in livestock and poultry genomes

[0052] 1. Select 5,000 Chinese Holstein cows born between 2018 and 2022 from multiple large ranches in Jiangsu, collect tail venous blood and extract genomic DNA. Use NanoDrop2000 to detect DNA purity, Qubit2.0 to accurately quantify DNA concentration, and agarose gel electrophoresis to detect genomic DNA integrity. The qualified DNA was sequenced for whole genome sequencing, and the sequencing result data set for each cow contained 14,264,452 SNP sites;

[0053] 2. Divide the haploid domain as the similarity interval based on the LD threshold:

[0054] (1) Calculate the LD value r between any two SNPs in the above 5000 cattle genome sequencing result data set according to formula (i): 2 :

[0055]

[0056] Where P(AB) is the frequency of AB gametes; P(A) and P(a) are the frequencies of alleles A and a respectively; P(B) and P(b) are the frequencies of alleles B and b respectively;

[0057] (2) The calculation results are arranged in order according to the position of the SNP sites in the genome to form a symmetric matrix H;

[0058] (3) Set the threshold R 1 =0.8, taking the first SNP site of chromosome 1 as the starting point, merging sites to the right, and calculating the minimum value of all elements in the square matrix formed by the starting SNP and the newly merged SNP each time in the matrix H, until the minimum value of all elements in the square matrix formed by the starting SNP and the newly merged SNP is less than the threshold;

[0059] (4) Remove the most recently merged SNP site, and use the starting SNP site and the last merged SNP site as the interval. If there are two or more SNP sites in the interval, the interval is regarded as a haploid region;

[0060] (5) Using the SNP closest to the right end of the newly divided haplotype block as a new starting point, repeat steps (3-4) until the genome is completely divided.

[0061] 3. Use the optimized K-means model to divide the SNP sites that are not included in the monomer domain interval to establish the similarity interval:

[0062] (1) The SNP sites not included in the monomer domain interval were integrated according to the position order and the genotypes were converted to 0, 1 or 2;

[0063] (2) The amount of genetic information of the above SNP site is calculated according to formula (ii), the conditional entropy of the above SNP site under the random variable site is calculated according to formula (iii), and the normalized information gain is calculated according to formula (iv):

[0064]

[0065] Among them, H(SNP x )、H(SNP x |SNP y ) and SU(SNP x ,SNP y ) are SNPs x Genetic information, SNP x In the random variable SNP y Conditional entropy and normalized information gain under P(SNP x =i) is SNP x The probability of genotype i;

[0066] (3) Using the modified Euclidean distance metric formula (v), calculate the Euclidean distance of the SNP sites not included in the monomer domain interval. i , calculate the site density value according to formula (vi):

[0067]

[0068] Among them, SNP xt and SNP yt is the genotype value of the corresponding t-th individual; SU(SNP x ,SNP y ) is the normalized information gain; e is the natural logarithm; r 1 is the effective density radius, which is usually replaced by the field radius; T is the total number of data points within the radius;

[0069] (4) According to the elbow rule, determine the K value at the inflection point where the sum of squared errors of all data suddenly changes with the number of clusters, and use this value as the optimal number of clusters K. 1 ;

[0070] (5) Select the top K with the largest density value 1 SNP sites were used as the initial cluster centers, and the distances between all SNP sites not included in the monomer domain interval and the cluster center vectors were calculated based on the improved Euclidean distance measurement formula. Based on this, each SNP site was classified into the category of the cluster center closest to it.

[0071] (6) Calculate K 1The sum of squared errors of the Euclidean distances between SNP sites in each cluster is calculated, and the SNP site with the smallest absolute value of the sum of squared errors of all SNP sites in each cluster is used as the new cluster center, and the distances between all SNPs in the genome and the new K are recalculated. 1 The Euclidean distances of cluster centers are re-divided into K 2 Clusters. Repeat this step until the cluster center no longer changes.

[0072] 4. Screening of high genetic information loci through variance analysis technology:

[0073] (1) For each divided interval in steps 2 and 3, convert the SNP site data into 0, 1 or 2 format, and calculate the genotype scaling value of each SNP site according to formula (vii), and the variance of each SNP value after scaling is equal to 1 and the mean value is equal to 0:

[0074]

[0075] in, and θ(SNP α ) are SNPs α The mean and standard deviation of the genotype values ​​of the locus in this Holstein cow population;

[0076] (2) Based on the genotype scaling value, the covariance matrix V between the SNP sites in each divided interval is calculated, and the SNP after removing the interference of all other SNP sites in the interval is calculated according to formula (viii). α Corrected value of genomic information explained by the site:

[0077]

[0078] Among them, V δδ 、V αα and V δα SNP δ and SNP α The variance and covariance between them, N is the number of SNP sites in the interval;

[0079] (3) Select S within each genome partition interval α The first three largest SNP sites are used as candidate high genetic information sites in the interval, and the (co)variance correction matrix between the remaining SNP sites in the interval is calculated according to formula (ix):

[0080]

[0081] Among them, V βγ Select SNPs for this interval ξ The scaling matrix before the site; V βξ 、Vγξ and V ξξ Other SNP sites and SNP in the interval are ξ The (co)variance size; N is the number of SNP sites in the interval. After correction, the next round of candidate high genetic information SNP sites in the interval are screened;

[0082] (4) After each round of screening of candidate high genetic information SNP sites, the total proportion of genomic information in each genome division interval is calculated according to formula (x):

[0083]

[0084] Among them, θ 0 2 is the total variance of the matrix before selecting the SNP site, θ t 2 is the residual variance of the matrix after selecting SNP sites, and t is the number of selected sites. For each divided interval, a threshold w=0.9 is set. When the total proportion of genomic information explained by the SNP sites selected in this interval reaches the threshold, the process stops and slides to the next genomic division interval for screening.

[0085] 5. After the 14,264,452 SNP sites contained in the sequencing result data set were screened, a set of 140K high genetic information SNP sites was obtained and output;

[0086] 6. For the 140K high genetic information SNP locus set, continue the above steps 1-5 to obtain a 10K high genetic information SNP locus set.

[0087] Example 2: Verification of the accuracy of interpolation of a set of highly informative SNP sites

[0088] 1. According to the sequencing result data set of each cow in Example 1, randomly select 140K site and 10K site to compare with the selected Holstein cattle 140K and 10K high genetic information SNP sites, respectively, and select the 100K chip site commonly used on the market to compare with the selected Holstein cattle 140K;

[0089] 2. Beagle (v4.0) was used to interpolate the above chip data, and Plink (v1.90) was used to perform quality control to eliminate sites with MAF < 0.05. The interpolation results are shown in Table 1. GraphPad Prism (v9.0) was used to make a bar graph of the chip interpolation results selected by different selection methods to more clearly see the difference in interpolation accuracy: Figure 2 As shown, the 10K sites selected by the new SNP selection method are compared with the 10K sites selected randomly; Figure 3As shown, the 140K loci selected by the new SNP selection method are compared with the 140K loci selected randomly and the 100K chip loci commonly used on the market.

[0090] 3. The interpolation accuracy of the selected 10K sites was good, with the overall level reaching 0.89±0.04 ( Figure 2 ), and its accuracy was higher than that of the randomly selected 10K loci on all 29 chromosomes. The 140K loci selected also had excellent interpolation accuracy, with an overall level of 0.86±0.02 ( Figure 3 ), which is significantly higher than the 0.83±0.02 of the commonly used 100K chip sites on the market, and also significantly higher than the 0.76±0.06 of the randomly selected 140K sites. The highest interpolation accuracy was for chromosome 19, with a value of 0.91, and on 29 chromosomes, the interpolation accuracy of the new SNP selection method was higher than that of the other two methods ( Figure 3 ).

[0091] Table 1 Statistical results of interpolation accuracy

[0092]

[0093] Embodiment 3: Device, computer-readable storage medium, and computer program product capable of implementing a method for screening a set of high genetic information SNP sites

[0094] Figure 4 An example of a physical structure diagram of an electronic device is shown in FIG. Figure 4 As shown, the electronic device may include: a processor (Processor) 410, a memory (Memory) 420, a communication interface (Communications Interface) 430 and a communication bus 440, wherein the processor 410, the memory 420, and the communication interface 430 communicate with each other through the communication bus 440. The processor 410 may call the logic instructions in the memory 420 to execute the high genetic information SNP site set screening method.

[0095] In addition, the logic instructions in the above-mentioned memory 420 can be implemented in the form of software functional units and when sold or used as independent products, can be stored in a computer-readable storage medium, including: USB flash drive, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk and other media that can store program codes.

[0096] On the other hand, the present invention also provides a computer program product, which includes a computer program. When executed by a processor, the computer can execute the above-mentioned method for screening a set of high genetic information SNP sites.

Claims

1. A method for screening a set of high genetic information SNP sites, characterized in that: The screening method comprises: (1) Divide haploid domains based on LD threshold; (2) The SNP sites not included in the monomer domain interval were divided by K-means model; (3) Using variance analysis technology, representative SNP sites within the similar genomic intervals divided in the previous step were screened out; (4) After the screening is completed, a set of SNP sites with high genetic information is output.

2. The screening method according to claim 1, characterized in that In step 1, the haploid domain is divided based on the LD threshold, including: (11) Calculate the LD value between any two SNP sites in the genome of each individual; (12) The calculation results are arranged in order according to the position of the SNP sites in the genome to form a symmetric matrix; (13) Setting the LD threshold, taking the first SNP site of the first chromosome as the starting point, merging sites to the right, and calculating the minimum value of all elements in the square matrix formed by each starting SNP and the newly merged SNP in the aforementioned matrix until the minimum value of all elements in the square matrix formed by the starting SNP and the newly merged SNP is less than the threshold; (14) Remove the most recently merged SNP site, and use the starting SNP site and the last merged SNP site as the interval. If there are two or more SNP sites in the interval, the interval is regarded as a haploid region; (15) Using the SNP closest to the right end of the newly divided haplotype block as a new starting point, repeat steps 13-14 until the genome is completely divided.

3. The screening method according to claim 2, characterized in that: The LD threshold in step 13 is a specific value between 0 and 1.

4. The screening method according to claim 1, characterized in that: In step 2, the SNP sites that are not included in the monomer domain interval are divided by the optimized K-means model, including: (21) SNP sites that were not included in the monomer domain interval were integrated in order of position; (22) The amount of genetic information of the above SNP site is calculated by formula 221, the conditional entropy under the random variable site is calculated by formula 222, and the normalized information gain is calculated by formula 223: Among them, H(SNP x )、H(SNP x |SNP y ) and SU(SNP x ,SNP y ) are SNPs x Genetic information, SNP x In the random variable SNP y Conditional entropy and normalized information gain under P(SNP x =i) is SNP x The probability of genotype i; (23) The Euclidean distance between the above SNP sites is calculated using the modified Euclidean distance metric formula 231, and the density distance value of all SNP sites is calculated using formula 232: Among them, SNP xt and SNP yt is the genotype value of the corresponding t-th individual; SU(SNP x ,SNP y ) is the normalized information gain; e is the natural logarithm; r1 is the effective density radius, which is generally replaced by the field radius; T is the total number of data points within the radius; (24) According to the elbow rule, the optimal number of clusters K1 is determined, and the first K1 SNP sites with the largest density values ​​calculated in step 23 are selected as the initial cluster centers. Based on the improved Euclidean distance measurement formula, the distances between all SNP sites not included in the monomer domain interval and the vectors of each cluster center are calculated again, and the SNP sites are classified into the category of the closest cluster center; (25) Calculate the sum of squared errors of the Euclidean distances between the SNP sites in the K1 clusters, take the SNP site with the smallest absolute value of the sum of squared errors among all SNP sites in each cluster as the new cluster center, recalculate the Euclidean distances between all SNP sites in the genome and the new K1 cluster centers and re-divide them into K2 clusters. Repeat this step until the cluster center no longer changes.

5. The screening method according to claim 1, characterized in that In step 3, the sites are screened by variance analysis techniques, including: (31) For each divided interval in step 1-2, the SNP site data is converted into 0, 1 or 2, and the genotype scaling value of each SNP site is calculated by formula 311: in, and θ(SNP α ) are SNPs α The mean and standard deviation of the genotype values ​​of the locus in the population; (32) Based on the genotype scaling value, the covariance matrix between the SNP sites in each divided interval is calculated, and the correction value of the amount of genomic information explained by each SNP site after removing the interference of other SNP sites in the interval is calculated by formula 321: Among them, V δδ 、V αα and V δα SNP δ and SNP α The variance and covariance between them, N is the number of SNP sites in the interval (33) Select the SNP site with the largest correction value of the amount of genome information explained in each genome division interval as the representative site of the interval, and use formula 331 to calculate the (co)variance correction matrix between the remaining SNP sites in the interval. After correction, screen the next round of candidate high genetic information SNP sites in the interval; Among them, V βγ Select SNPs for this interval ξ The scaling matrix before the site; V βξ 、V γξ and V ξξ Other SNP sites and SNP in the interval are ξ The (co)variance size, N is the number of SNP sites in this interval; (34) After each round of screening of candidate high genetic information SNP sites, the total proportion of genomic information in each genome division interval is calculated by formula 341: Among them, θ0 2 is the total variance of the matrix before selecting the SNP site, θ t 2 is the residual variance of the matrix after selecting SNP sites, t is the number of selected sites, and for each divided interval, a threshold w is set. When the total proportion of genomic information explained by the SNP sites selected in this interval reaches the threshold, it stops and slides to the next genome division interval for screening.

6. The screening method according to claim 5, characterized in that In step 32, one or more SNP sites with the largest genomic information correction value in each genomic division interval are selected.

7. The screening method according to claim 5, characterized in that The threshold w in step 34 is a specific value between 0.6 and 1.

8. A computer device comprising a memory, a processor and a computer program stored in the memory, characterized in that: When the processor executes the computer program, the screening method according to any one of claims 1 to 7 is implemented.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the screening method according to any one of claims 1 to 7 is implemented.

10. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the screening method according to any one of claims 1 to 7 is implemented.

Images