A snowflake black cattle breeding method and system based on whole genome association analysis

Abstract

The application provides a snowflake black cattle breeding method and system based on whole genome association analysis, relates to the technical field of black cattle breeding, and comprises the following steps: extracting PacBio long read HiFi sequencing data and second-generation whole genome resequencing data in a blood sample, and extracting Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data in a tissue sample; performing genome assembly to obtain a chromosome-level reference genome sequence; performing gene structure annotation on the reference genome sequence; detecting a whole genome SNP marker data set of a breeding population; screening whole genome SNP marker data meeting a preset significance threshold; performing functional annotation to obtain molecular markers for snowflake black cattle breeding. The application provides reliable and practically functional molecular marker support for the directional breeding of snowflake black cattle meat quality-related key economic traits.

Description

Technical Field

[0001] This invention relates to the field of black cattle breeding technology, and more specifically, to a method and system for breeding Snowflake Black Cattle based on genome-wide association analysis. Background Technology

[0002] Snowflake beef, also known as marbled beef, gets its name from the white fat evenly distributed between red muscle fibers, resembling snowflakes (with an intramuscular fat content of over 13%). It not only boasts a superior taste and rich flavor but is also rich in conjugated linoleic acid (CLA) and beneficial components such as calcium, iron, potassium, zinc, phosphorus, and various vitamins and minerals. The meat is tender, juicy, and extremely nutritious, making it highly sought after both domestically and internationally, commanding a high price. However, in China, due to the low number of purebred Wagyu cattle and slow breeding rates, a significant industry has not yet been established. Consequently, domestic production of snowflake beef relies primarily on crossbreeding, using purebred Wagyu as the sire and other breeds as the dam to breed F1 hybrid Wagyu. Statistics show that only 30%-40% of F1 fattened and slaughtered beef reaches A3 grade or higher, while the proportion of purebred Wagyu fattened and slaughtered reaches over 85%. Each additional grade of beef is expected to increase the price by 30%, highlighting the significant economic difference between crossbred and purebred Wagyu cattle. In addition to selecting appropriate breeds, it is also necessary to establish a sound breeding system for producing Wagyu beef.

[0003] The current Wagyu cattle population in China suffers from low breeding efficiency, a lack of core breeds, uneven production efficiency, and slow genetic progress, severely hindering the development of the high-end Wagyu beef industry. Genome-wide selection (WHM) technology is a current hot topic in livestock breeding research, allowing for accurate early selection of individuals without relying on phenotypic information, thereby significantly shortening generation intervals and accelerating population genetic progress. It is currently widely used in dairy cattle breeding. However, due to the complexity of genomic regulation during Wagyu growth and development, existing Wagyu genome sequences contain many gaps, particularly complex repetitive sequences, making in-depth genomic feature analysis difficult. Furthermore, accurately and effectively mining the hidden genetic information within the Wagyu genome remains a common challenge for the industry, hindering the application of WHM in Wagyu breeding.

[0004] No effective solutions have yet been proposed to address the problems in the relevant technologies. Summary of the Invention

[0005] In view of this, the present invention provides a method and system for breeding Snowflake Black Cattle based on genome-wide association analysis to solve the aforementioned problems.

[0006] To solve the above problems, the specific technical solution adopted by the present invention is as follows:

[0007] According to one aspect of the present invention, a method for breeding Snowflake Black Cattle based on genome-wide association analysis is provided, the method comprising the following steps:

[0008] S1. Collect blood and tissue samples from the Snowflake Black Cattle breeding population, and extract PacBio long-read HiFi sequencing data and second-generation whole-genome resequencing data from the blood samples, as well as Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data from the tissue samples.

[0009] S2. Based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data, genome assembly was performed to obtain a chromosome-level reference genome sequence; based on multi-tissue full-length transcriptome sequencing data, gene structure annotation was performed on the reference genome sequence to obtain gene set and transcript information;

[0010] S3. Align the second-generation whole-genome resequencing data to the reference genome sequence and detect the whole-genome SNP marker dataset of the breeding population;

[0011] S4. Collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and combine them with the whole genome SNP marker dataset. Through whole genome association analysis, screen whole genome SNP marker data that meet the preset significance threshold.

[0012] S5. Using gene set and transcript information, functional annotation was performed on the screened whole-genome SNP marker data to obtain molecular markers for the breeding of Snowflake Black Cattle.

[0013] Preferably, the process of collecting standardized phenotypic data from all individuals in the Snowflake Black Cattle breeding population, and combining this data with a genome-wide SNP marker dataset, to screen genome-wide SNP marker data that meet a preset significance threshold through genome-wide association analysis includes the following steps:

[0014] S41. Collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and obtain corrected unbiased phenotypic values ​​by correcting the standardized phenotypic data.

[0015] S42. Based on the whole genome SNP marker dataset, construct an SNP relationship network, and use local co-optimization to divide the SNP relationship network to obtain several SNP functional modules.

[0016] S43. For each SNP functional module, the association significance between the joint effect of all SNPs and the corrected unbiased phenotypic value is calculated using the multi-label association analysis method to obtain the association value of each SNP functional module, and the SNP functional modules that are significantly associated are selected according to the preset significance threshold.

[0017] S44. For each significantly associated SNP functional module, and using the transformation of the cumulative distribution function, select several SNPs as representative markers in each SNP functional module, and construct the screened whole genome SNP marker data from all representative markers.

[0018] Preferably, the step of constructing an SNP relationship network based on a whole-genome SNP marker dataset and dividing the SNP relationship network into several SNP functional modules using a local co-optimization approach includes the following steps:

[0019] S421. Based on the whole genome SNP marker dataset, calculate the biological functional association coefficient between any two SNP marker datasets;

[0020] S422. Construct an SNP relationship network with SNP-labeled data as nodes and biological function correlation coefficients as edge weights;

[0021] S423. By calculating the validity of each node in the SNP relationship network and combining it with local collaborative optimization, the SNP relationship network is divided into several SNP functional modules.

[0022] Preferably, the calculation of the biological functional association coefficient between any two SNP marker data based on the whole genome SNP marker dataset includes the following steps:

[0023] S4211. Based on the Hi-C chromosome conformation, a chromosome-level chromatin interaction matrix is ​​constructed using the captured sequencing data, and the normalized contact frequency between any two SNP marker data regions is extracted as the spatial interaction strength coefficient.

[0024] S4212. Based on gene set and transcript information, combined with a pre-set biological pathway database, pathway annotation is performed on each SNP marker data, and semantic similarity is used to calculate the functional similarity coefficient between any two SNP marker data.

[0025] S4213. Perform quantitative trait locus localization analysis on multi-tissue full-length transcriptome sequencing data to obtain the target genes associated with each SNP as a regulatory site and their effect direction and intensity, and calculate the correlation coefficient of any two SNP marker data on the regulatory expression profile based on the effect direction and intensity.

[0026] S4214. The spatial interaction strength coefficient, functional similarity coefficient, and correlation coefficient are weighted and summed according to preset weight coefficients to obtain the biological functional correlation coefficient.

[0027] Preferably, the step of dividing the SNP relationship network into several SNP functional modules by calculating the effectiveness of each node in the SNP relationship network and combining it with local collaborative optimization includes the following steps:

[0028] S4231. Mark all nodes as unvisited and initialize the functional module list to empty;

[0029] S4232. Calculate the validity of each node in the SNP relationship network based on the sum of the biological functional correlation coefficients of all edges connected to the node, and select the node with the highest validity as the core node.

[0030] S4233. Starting from the core node, perform local collaborative expansion processing to generate candidate SNP functional modules, and mark the core node as visited.

[0031] S4234. Repeat steps S4232 to S4233 until all nodes are marked as visited, and obtain the preliminary SNP functional module partitioning results.

[0032] S4235. Analyze the conflict nodes in the preliminary SNP functional module partitioning results, and based on the analysis results, reallocate the preliminary SNP functional module partitioning results to obtain the final SNP functional module partitioning results.

[0033] Preferably, the step of performing local collaborative expansion processing starting from the core node to generate candidate SNP functional modules and marking the core node as visited includes the following steps:

[0034] S42331. Based on the biological functional correlation coefficient between the core node and its neighboring nodes, calculate the correlation strength between the core node and all its neighboring nodes, and sort the neighboring nodes in descending order of correlation strength.

[0035] S42332. Initialize the candidate SNP function module, mark the core node as visited and include it in the candidate module, select the first-ranked neighbor node to add to the candidate SNP function module and mark it as visited, integrate the neighbor nodes of all nodes in the candidate SNP function module to form a set of neighbors to be expanded, and remove visited nodes at the same time.

[0036] S42333. Calculate the internal and external degrees of the current candidate SNP functional modules, and calculate the fitness of each node in the neighbor set to be expanded based on the internal and external degrees.

[0037] S42334. Add neighbor nodes with fitness greater than zero to the candidate SNP function module and mark them as visited. If no node in the neighbor set to be expanded satisfies fitness greater than zero, end the local collaborative expansion process.

[0038] Preferably, for each SNP functional module, the step of calculating the significance of the association between the joint effect of all SNPs and the corrected unbiased phenotypic value using a multi-label association analysis method to obtain the association value of each SNP functional module, and selecting the SNP functional modules with significant association according to a preset significance threshold, includes the following steps:

[0039] S431. Extract the genotype matrix of all SNP marker data in each SNP functional module and match it with the corrected unbiased phenotypic values ​​to construct the association analysis dataset for each SNP functional module.

[0040] S432. Calculate the significance of the association between the joint effect of all SNPs within each SNP functional module and the unbiased phenotypic value using the multi-label association analysis method, and obtain the association value of each SNP functional module.

[0041] S433. Compare the correlation values ​​of SNP functional modules with a preset significance threshold, and select SNP functional modules with correlation values ​​less than the significance threshold as significantly correlated SNP functional modules.

[0042] Preferably, the step of extracting the genotype matrix of all SNP marker data within each SNP functional module and matching it with the corrected unbiased phenotypic values ​​to construct the association analysis dataset for each SNP functional module includes the following steps:

[0043] S4311. For each SNP functional module, extract the genotype information of all SNP marker data in each individual of the breeding population, and convert the genotypes into numerical data using numerical encoding to construct an initial genotype matrix.

[0044] S4312. Perform quality control on the initial genotype matrix by removing SNP markers with genotype deletion rates exceeding a preset threshold and filling the remaining missing values ​​with the mean to obtain a complete genotype matrix.

[0045] S4313. Match the complete genotype matrix with the corrected unbiased phenotypic values ​​to form the association analysis dataset for each SNP functional module.

[0046] Preferably, the step of calculating the significance of the association between the joint effect of all SNPs within each SNP functional module and the unbiased phenotypic value using multi-label association analysis to obtain the association value for each SNP functional module includes the following steps:

[0047] S4321. Based on the association analysis dataset of each SNP functional module, select a multi-marker statistical model suitable for association analysis of high-dimensional genetic data from the preset multi-marker statistical model library.

[0048] S4322. Based on the selected multi-label statistical model, calculate the joint effect statistic between all SNP labels and unbiased phenotypic values ​​in the SNP functional module and their corresponding original association values.

[0049] S4323. Perform multiple hypothesis testing to correct the original association values, and use the corrected original association values ​​as the final association values ​​for each SNP functional module.

[0050] According to another aspect of the present invention, a Snowflake Black Cattle breeding system based on genome-wide association analysis is provided, the system comprising:

[0051] The data acquisition module is used to collect blood and tissue samples from the Snowflake Black Cattle breeding population, and extract PacBio long-read HiFi sequencing data and second-generation whole-genome resequencing data from the blood samples, as well as extract Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data from the tissue samples.

[0052] The assembly sequencing module is used to assemble the genome based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data to obtain a chromosome-level reference genome sequence; and to annotate the reference genome sequence with gene structure based on multi-tissue full-length transcriptome sequencing data to obtain gene set and transcript information.

[0053] The data alignment module is used to align second-generation whole-genome resequencing data to a reference genome sequence and to detect whole-genome SNP marker datasets of the breeding population.

[0054] The data filtering module is used to collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and combined with the whole genome SNP marker dataset, to filter whole genome SNP marker data that meet the preset significance threshold through genome-wide association analysis;

[0055] The functional annotation module is used to perform functional annotation on the screened whole-genome SNP marker data using gene set and transcript information, so as to obtain molecular markers for the breeding of Snowflake Black Cattle.

[0056] The beneficial effects of this invention are as follows:

[0057] 1. This invention assembles a chromosome-level reference genome specific to Snowflake Black Cattle and completes high-precision gene structure annotation. Based on this specific reference genome, it enables accurate detection of SNP markers in the entire genome of the breeding population. Combined with gene set and transcript information, it completes functional annotation, breaking through the problems of SNP detection bias and inaccurate annotation caused by traditional breeding that relies on general bovine reference genomes. From genome assembly, marker detection to association screening and functional annotation, it forms a complete process of specific breeding system adapted to Snowflake Black Cattle, which greatly improves the accuracy, scientificity and targeting of Snowflake Black Cattle molecular breeding. It provides reliable and practically functional molecular marker support for the targeted breeding of key economic traits related to Snowflake Black Cattle beef quality.

[0058] 2. This invention integrates multi-dimensional information such as three-dimensional genome spatial interactions, similar biological pathway functions, and gene expression regulation associations to calculate the biological functional association coefficients between SNPs, thereby constructing an SNP relationship network that closely reflects the essence of biological functions. Furthermore, it utilizes a local collaborative optimization method to select core nodes based on effectiveness and to orderly expand and divide functionally homogeneous SNP modules according to association strength and fitness, overcoming the limitations of traditional genetic linkage-based division. Multiple hypothesis testing is used to reduce false positives, and after screening for significantly associated modules, representative markers are selected using the cumulative distribution function, significantly improving the accuracy and reliability of association analysis results. The final selected SNP markers possess both genetic significance and biological functional effectiveness, making marker screening for the molecular breeding of Wagyu Black Cattle more targeted and scientific, providing precise and functionally valuable molecular marker support for the subsequent efficient breeding of high-quality Wagyu Black Cattle. Attached Figure Description

[0059] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly described below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings:

[0060] Figure 1 This is a flowchart of a Snowflake Black Cattle breeding method based on genome-wide association analysis according to an embodiment of the present invention;

[0061] Figure 2 This is a schematic diagram of a Snowflake Black Cattle breeding system based on genome-wide association analysis according to an embodiment of the present invention.

[0062] In the picture:

[0063] 1. Data acquisition module; 2. Assembly and sequencing module; 3. Data alignment module; 4. Data filtering module; 5. Functional annotation module. Detailed Implementation

[0064] To enable those skilled in the art to better understand the technical solutions in this application, the technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this application.

[0065] According to an embodiment of the present invention, a method and system for breeding Snowflake Black Cattle based on genome-wide association analysis is provided.

[0066] The present invention will now be further described in conjunction with the accompanying drawings and specific embodiments, such as... Figure 1 As shown, according to an embodiment of the present invention, a method for breeding Snowflake Black Cattle based on genome-wide association analysis is provided, the method comprising the following steps:

[0067] S1. Collect blood and tissue samples from the Snowflake Black Cattle breeding population, and extract PacBio long-read HiFi sequencing data and second-generation whole-genome resequencing data from the blood samples, as well as Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data from the tissue samples.

[0068] It should be noted that the Snowflake Black Cattle breeding population consisted of healthy, representative individuals. Blood samples, as well as samples from multiple tissues including muscle, fat, and liver, were collected. Blood samples were used to extract high-molecular-weight genomic DNA, while tissue samples were rapidly cryopreserved in liquid nitrogen or RNA preservation solution immediately after being removed from the body.

[0069] In the acquisition of PacBio long-read HiFi sequencing data, high molecular weight genomic DNA was extracted from blood samples using the phenol-chloroform method or a commercially available large-fragment DNA extraction kit. The purity and integrity of the DNA were detected by pulsed-field electrophoresis, requiring no significant degradation of the DNA, an optical density 260 / 280 ratio between 1.8 and 2.0, and a fragment length ≥50 kbps. A high-fidelity library was constructed, and the DNA was fragmented, damaged, and end-repaired before being ligated to SMRTbell adapters. Exonuclease digestion was then performed to remove unligated fragments. On the sequencing platform, a high-fidelity sequencing mode was used, with a sequencing depth of no less than 30× per SMRT cell, to obtain long-read high-fidelity reads with high accuracy (>99.9%).

[0070] Second-generation whole-genome resequencing data was also obtained by extracting genomic deoxyribonucleic acid (DNA) from blood samples, requiring a total DNA amount ≥1 μg and a concentration ≥20 nanograms / µL, and constructing a paired-end sequencing library with 350-450 base pairs of insert fragments. Paired-end 150-base-pair sequencing was performed on the DNBSEQ-T7 platform to ensure uniform coverage of the entire genome and accuracy of single nucleotide polymorphism (SNP) detection.

[0071] The acquisition of Hi-C chromosome conformation capture sequencing data involves taking frozen tissue samples, grinding them into powder in liquid nitrogen, and then cross-linking them with formaldehyde to fix the spatial conformation of chromatin. The nucleus is then released by adding lysis buffer. Restriction endonucleases are used to digest the chromatin, producing sticky ends. These ends are filled with biotin-labeled nucleotides, and blunt-end ligation is performed to link spatially adjacent DNA fragments into chimeric molecules. After reversing the cross-linking, the DNA is purified to remove unlinked biotin, and the DNA fragments are fragmented to 300-500 base pairs. The biotin-labeled ligation products are enriched using streptavidin magnetic beads, and an Illumina sequencing library is constructed. Paired-end 150-base-pair sequencing is then performed on the Illumina NovaSeq platform.

[0072] Acquiring full-length transcriptome sequencing data from multiple tissues, including muscle, fat, and liver, involved collecting frozen samples and extracting total ribonucleic acid (RNA) using a RNA extraction kit. RNA integrity was assessed by agarose gel electrophoresis, requiring an RNA integrity index ≥8. Messenger RNA was purified from the total RNA and reverse transcribed to synthesize complementary deoxyribonucleic acid (DNA). After polymerase chain reaction (PCR) amplification, the DNA was fractionated by size to construct an SMRTbell library. Sequencing was then performed on the Revio platform to obtain full-length reads covering the complete transcriptome.

[0073] The key parameter settings for PacBio long-read HiFi sequencing data, second-generation whole-genome resequencing data, Hi-C chromosome conformation capture sequencing data, and multi-tissue full-length transcriptome sequencing data are shown in Table 1.

[0074] Table 1 Key parameters of sequencing data

[0075] Data types sequencing platform sequencing mode Sequencing depth Key parameters Data quality requirements PacBio long-read HiFi sequencing PacBio Sequel II HiFi mode ≥30× Read length ≥ 10kb, error rate < 1% Genome assembly error rate <1.5% Second-generation whole-genome resequencing DNBSEQ-T7 150bp dual-ended ≥30× Insert fragment 350-450bp SNP detection accuracy ≥ 95% Hi-C chromosome conformation capture sequencing Illumina NovaSeq 150bp dual-ended ≥50× 100kb resolution Contact matrix sparsity <15% Multi-tissue full-length transcriptome sequencing Revio Full-length transcript ≥50× RIN≥8.0 Transcript integrity ≥95%

[0076] S2. Based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data, genome assembly was performed to obtain a chromosome-level reference genome sequence; based on multi-tissue full-length transcriptome sequencing data, gene structure annotation was performed on the reference genome sequence to obtain gene set and transcript information;

[0077] It should be noted that the genome assembly based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data to obtain a chromosome-level reference genome sequence includes: quality filtering of the PacBio long-read HiFi sequencing data to remove low-quality reads shorter than 5kb, obtaining high-accuracy, high-fidelity circularized consistent sequencing reads; then, using genome assembly software, de novo assembly of the filtered HiFi reads to construct genome contigs, and preliminary evaluation of the assembly results to ensure the continuity and accuracy of the contigs; aligning the Hi-C chromosome conformation capture sequencing data to the preliminarily assembled contig sequences to identify chromatin spatial interactions between different contigs; using the interaction frequency matrix of the Hi-C data, using chromosome mounting tools to cluster, sort, and orient the contigs, constructing a chromosome-level scaffold sequence; finally, evaluating the integrity of the chromosome-level assembly results and verifying the mounting accuracy using a Hi-C interaction heatmap, ultimately obtaining a chromosome-level reference genome sequence.

[0078] The process involves using multi-tissue full-length transcriptome sequencing data to annotate the reference genome sequence, obtaining gene sets and transcript information. This includes: initially aligning the multi-tissue full-length transcriptome sequencing data to the chromosome-level reference genome sequence to obtain the precise location information of each transcript in the genome, including exon regions and splicing sites; performing cluster analysis on the initially aligned transcripts, merging identical or highly similar transcripts based on genomic coordinates and splicing patterns, removing redundant transcripts, and obtaining a non-redundant full-length transcript set; and based on the non-redundant transcript set, combining it with gene models of known reference species. The model is developed by correcting the gene structure on the reference genome, locating the transcription start site, exon-intron boundary, and transcription termination site. New transcripts and alternative splicing events inconsistent with known gene models are then identified to supplement and improve the gene structure information of the reference genome, resulting in a corrected gene model. Finally, coding region prediction and protein function annotation are performed on the corrected gene model, and the biological function of the genes is inferred by comparing it with functional databases. Ultimately, a gene set containing gene location, structural information, and functional descriptions, as well as transcript information including all non-redundant transcript sequences and alternative splicing patterns, are constructed.

[0079] S3. Align the second-generation whole-genome resequencing data to the reference genome sequence and detect the whole-genome SNP marker dataset of the breeding population;

[0080] It should be noted that when aligning the second-generation whole-genome resequencing data to the reference genome sequence, high-performance alignment software, such as BWA-MEM or Bowtie2, is used to accurately align the quality-controlled second-generation whole-genome resequencing reads with the constructed chromosome-level reference genome sequence. The alignment results are sorted and deduplicated, and population-level detection of single nucleotide polymorphisms (SNPs) is performed using variation detection tools such as GATK or SAMtools. By setting filtering criteria such as base quality value ≥20, sequencing depth ≥4×, and minor allele frequency ≥0.05, a high-confidence SNP marker dataset covering the entire genome is finally obtained.

[0081] S4. Collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and combine them with the whole genome SNP marker dataset. Through whole genome association analysis, screen whole genome SNP marker data that meet the preset significance threshold.

[0082] In a preferred embodiment, the collection of standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, combined with a genome-wide SNP marker dataset, and the screening of genome-wide SNP marker data that meet a preset significance threshold through genome-wide association analysis includes the following steps:

[0083] S41. Collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and obtain corrected unbiased phenotypic values ​​by correcting the standardized phenotypic data.

[0084] It should be noted that meat quality phenotypic data of all individuals in the Wagyu Black Cattle breeding population were collected, including but not limited to key economic traits such as marbling grade, tenderness, backfat thickness, and eye muscle area. After outlier removal and normality testing of the meat quality phenotypic data, the phenotypic data were corrected. Non-genetic factors such as environmental effects, feeding batch effects, sex effects, and age effects were treated as fixed effects, and the kinship matrix between individuals was treated as a random effect for regression analysis. The influence of these non-genetic factors on the phenotypic values ​​was eliminated, and the residuals were extracted as the corrected unbiased phenotypic values ​​for subsequent genome-wide association analysis, thereby eliminating environmental interference and improving the accuracy of genetic effect assessment.

[0085] S42. Based on the whole genome SNP marker dataset, construct an SNP relationship network, and use local co-optimization to divide the SNP relationship network to obtain several SNP functional modules.

[0086] As a preferred embodiment, the step of constructing an SNP relationship network based on a whole-genome SNP marker dataset and dividing the SNP relationship network into several SNP functional modules using a local co-optimization approach includes the following steps:

[0087] S421. Based on the whole genome SNP marker dataset, calculate the biological functional association coefficient between any two SNP marker datasets;

[0088] In a preferred embodiment, calculating the biological functional association coefficient between any two SNP marker datasets based on a whole-genome SNP marker dataset includes the following steps:

[0089] S4211. Based on the Hi-C chromosome conformation, a chromosome-level chromatin interaction matrix is ​​constructed using the captured sequencing data, and the normalized contact frequency between any two SNP marker data regions is extracted as the spatial interaction strength coefficient.

[0090] It should be noted that the obtained Hi-C sequencing data underwent quality assessment and filtering to remove low-quality reads and PCR repeats. Using dedicated Hi-C analysis software, the filtered Hi-C reads were aligned to the constructed chromosome-level reference genome sequence. Valid intrachromatin and interchromatin interaction pairs were identified based on the alignment positions. The genome was divided into fixed-size windows, such as 100kb or 500kb, and a genome-wide chromatin interaction matrix was constructed by counting the number of primitive interactions between pairs within each window.

[0091] Furthermore, to eliminate the influence of factors such as sequencing depth and genome accessibility, iterative correction or matrix balancing methods were used to normalize the interaction matrix, resulting in a standardized contact frequency matrix. For any two SNP markers, their physical locations on the reference genome were located, their respective windows were determined, and the standardized contact frequencies between corresponding windows were extracted from the standardized interaction matrix as the spatial interaction strength coefficients between the two SNPs.

[0092] S4212. Based on gene set and transcript information, combined with a pre-set biological pathway database, pathway annotation is performed on each SNP marker data, and semantic similarity is used to calculate the functional similarity coefficient between any two SNP marker data.

[0093] It should be noted that, based on the obtained gene set and transcript information, the genes located at or near each SNP marker are identified, such as those located within a gene, in a promoter region, or associated with an eQTL. These genes are then compared with a pre-defined biological pathway database to obtain the pathway set involved by each SNP-associated gene. For any two SNPs, pathway annotation sets of their associated genes are obtained, and semantic similarity is used to calculate the similarity between the two pathway sets. This semantic similarity is based on the pathway's hierarchical structure and information content, quantifying the degree of similarity between the two SNPs at the biological function level.

[0094] S4213. Perform quantitative trait locus localization analysis on multi-tissue full-length transcriptome sequencing data to obtain the target genes associated with each SNP as a regulatory site and their effect direction and intensity, and calculate the correlation coefficient of any two SNP marker data on the regulatory expression profile based on the effect direction and intensity.

[0095] It should be noted that quantitative analysis of full-length transcriptome sequencing data from multiple tissues can obtain the expression levels of each gene in each tissue sample, such as TPM or FPKM values. Combined with the detected whole-genome SNP marker dataset, eQTL localization analysis software is used to perform a genome-wide eQTL scan to analyze the significant association between each SNP and the expression levels of each gene. For each SNP, a set of target genes significantly associated with it (e.g., FDR < 0.05) is selected, and the effect direction and intensity of each target gene are recorded. For any two SNPs, their regulatory effect vectors on all target genes are constructed, with each vector element representing the effect value for each target gene. Genes without association are assigned a value of 0. The Pearson correlation coefficient or Spearman correlation coefficient between these two vectors is calculated as the correlation coefficient between the two SNPs on the regulatory expression profile.

[0096] S4214. The spatial interaction strength coefficient, functional similarity coefficient, and correlation coefficient are weighted and summed according to preset weight coefficients to obtain the biological functional correlation coefficient.

[0097] Specifically, since the three coefficients obtained have different dimensions and numerical ranges, they need to be normalized before weighted summation, mapping them to the 0-1 interval. Normalization methods can include min-max standardization or quantile standardization based on cumulative distribution. Based on the importance of the three dimensions in SNP functional association, the preset weight coefficients are: spatial interaction weight 0.3, functional similarity weight 0.3, and regulatory correlation weight 0.4, or objective weighting methods such as entropy weighting or principal component analysis can be used to determine the weights. The three normalized coefficients are multiplied by their corresponding weight coefficients and then summed to obtain a comprehensive biological functional association coefficient.

[0098] S422. Construct an SNP relationship network with SNP-labeled data as nodes and biological function correlation coefficients as edge weights;

[0099] It should be noted that each SNP in the whole-genome SNP marker dataset is treated as a node in the network. The biofunctional association coefficient between any two SNPs is used as the weight of the edge connecting these two nodes, constructing an undirected weighted network, namely the SNP relationship network. In this network, nodes represent SNP markers, and edges represent the functional associations between SNPs. The larger the edge weight, the stronger the functional association between the two SNPs at the spatial conformation, biological pathway, or transcriptional regulation level. The network is stored and represented in the form of an adjacency matrix, where the rows and columns of the matrix correspond to all SNP nodes, and the matrix elements are the biofunctional association coefficients of the corresponding SNP pairs.

[0100] S423. By calculating the validity of each node in the SNP relationship network and combining it with local collaborative optimization, the SNP relationship network is divided into several SNP functional modules.

[0101] As a preferred embodiment, the step of dividing the SNP relationship network into several SNP functional modules by calculating the effectiveness of each node in the SNP relationship network and combining it with local collaborative optimization includes the following steps:

[0102] S4231. Mark all nodes as unvisited and initialize the functional module list to empty;

[0103] It should be noted that before partitioning the SNP relationship network, the state of all nodes in the network needs to be initialized by marking all nodes as "unvisited," indicating that these nodes have not yet been assigned to any functional module. Simultaneously, an empty list of functional modules is initialized to store the various SNP functional modules obtained later.

[0104] S4232. Calculate the validity of each node in the SNP relationship network based on the sum of the biological functional correlation coefficients of all edges connected to the node, and select the node with the highest validity as the core node.

[0105] It's important to note that in a SNP network, the validity of each node is defined as the sum of the biofunctional association coefficients of all edges connected to that node, i.e., the node's weighted degree. Validity reflects the total strength of the SNP's functional associations within the overall network: a higher validity indicates a broader and closer functional association between the SNP and other SNPs, placing it in a more central position within the network. Therefore, the node with the highest validity from the set of unvisited nodes is selected as the core node, serving as the starting point for constructing a new SNP functional module.

[0106] S4233. Starting from the core node, perform local collaborative expansion processing to generate candidate SNP functional modules, and mark the core node as visited.

[0107] In a preferred embodiment, the step of performing local collaborative expansion processing starting from the core node to generate candidate SNP functional modules and marking the core node as visited includes the following steps:

[0108] S42331. Based on the biological functional correlation coefficient between the core node and its neighboring nodes, calculate the correlation strength between the core node and all its neighboring nodes, and sort the neighboring nodes in descending order of correlation strength.

[0109] Specifically, in the initial stage of local expansion, all neighboring nodes of the core node are obtained, i.e., the SNPs directly connected to the core node. For each neighboring node, the association strength with the core node is the biofunctional association coefficient between them, i.e., the edge weight. These neighboring nodes are then sorted in descending order of association strength to form a priority queue.

[0110] S42332. Initialize the candidate SNP function module, mark the core node as visited and include it in the candidate module, select the first-ranked neighbor node to add to the candidate SNP function module and mark it as visited, integrate the neighbor nodes of all nodes in the candidate SNP function module to form a set of neighbors to be expanded, and remove visited nodes at the same time.

[0111] It should be noted that initializing an empty candidate SNP functional module allows you to add core nodes to the module and mark them as visited.

[0112] The process involves selecting the top-ranked neighbor node from the sorted list, adding it to the candidate module, and marking it as visited. After adding these two nodes, the neighbor nodes of all nodes in the current candidate module are collected and merged into a set of neighbors to be expanded. When forming this set, nodes already marked as visited (i.e., core nodes and the added neighbor nodes themselves) need to be removed to ensure that subsequent nodes examined are all unvisited potential candidate nodes.

[0113] S42333. Calculate the internal and external degrees of the current candidate SNP functional modules, and calculate the fitness of each node in the neighbor set to be expanded based on the internal and external degrees.

[0114] It should be noted that, in order to determine whether each node in the neighbor set to be expanded is suitable for joining the current candidate SNP functional module, a fitness function is defined. This is achieved by calculating the internal degree of the current candidate SNP functional module. internal degree Defined as twice the sum of the biofunctional correlation coefficients of all edges connecting nodes within a candidate SNP functional module; simultaneously, the module's externality is calculated. , external degree Defined as the sum of the biofunctional correlation coefficients of all edges connecting nodes within a module to nodes outside the module.

[0115] Based on internality and externality, the compactness function of candidate SNP functional modules is defined as follows: ;

[0116] In the formula, α is a parameter adaptively adjusted based on the average clustering coefficient of the network, used to control the size of the module. For each node v in the neighbor set to be expanded, the change in its density before and after joining the current module is calculated: This change value represents the fitness of node v to the current module. A fitness greater than 0 indicates that adding node v can improve the overall density of the module and is suitable for inclusion; a fitness less than or equal to 0 indicates that adding node v will reduce the density of the module and is not suitable for inclusion.

[0117] S42334. Add neighbor nodes with fitness greater than zero to the candidate SNP function module and mark them as visited. If no node in the neighbor set to be expanded satisfies fitness greater than zero, end the local collaborative expansion process.

[0118] Specifically, based on the calculated fitness of each neighbor node, all nodes with a fitness greater than 0 are added to the candidate SNP functional module and marked as visited. The addition of these nodes increases the module's density and is a beneficial extension of the SNP functional module. After one round of additions, the set of neighbors to be expanded needs to be updated: all unvisited neighbors of the newly added nodes are added to the set, and the nodes are reordered according to the strength of their association with the core nodes, so that the nodes with the strongest associations are prioritized in the next iteration.

[0119] S4234. Repeat steps S4232 to S4233 until all nodes are marked as visited, and obtain the preliminary SNP functional module partitioning results.

[0120] S4235. Analyze the conflict nodes in the preliminary SNP functional module partitioning results, and based on the analysis results, reallocate the preliminary SNP functional module partitioning results to obtain the final SNP functional module partitioning results.

[0121] It should be noted that during the initial partitioning process, some nodes located at the boundaries of SNP functional modules may simultaneously have high functional association strength with the core or internal nodes of multiple SNP functional modules. This leads to them being repeatedly attempted to be included during the expansion of multiple modules; these nodes are called conflicting nodes. For these conflicting nodes, reallocation is required to optimize the final module partitioning. Specifically, for each conflicting node, its module affinity with each associated initial module is calculated. The module affinity calculation comprehensively considers the attraction of the core node of the SNP functional module to the conflicting node, as well as the attraction of all neighboring nodes connected to the node within the SNP functional module, and is obtained through weighted summation.

[0122] It should be understood that core node attractiveness refers to the attraction of the core node of the candidate SNP functional module to conflicting nodes. As the starting node in the formation of the SNP functional module, the core node represents the core functional direction of the SNP functional module, and the biofunctional correlation coefficient between the core node and the conflicting node directly reflects the degree of functional association between the two nodes. The greater the attractiveness of the core node to the conflicting node, the more likely the conflicting node is to be associated with the core function of the module.

[0123] Neighbor node attractiveness refers to the sum of the attractiveness of all neighboring nodes directly connected to the conflicting node within the candidate SNP functional module. The biofunctional association coefficient between each neighboring node and the conflicting node reflects the strength of the local functional association between the neighboring node and the conflicting node. Summing up the attractiveness of all neighboring nodes reflects the overall integration degree of the conflicting node with the internal structure of the entire SNP functional module. If a conflicting node has strong associations with multiple nodes within the module, it is more likely to belong to that SNP functional module.

[0124] In the specific calculation, different weight coefficients are assigned to the attractiveness of core nodes and the attractiveness of neighboring nodes within the SNP functional module, such as a weight of 0.4 for core nodes and 0.6 for neighboring nodes. These are then weighted and summed to obtain the final module affinity. The weight coefficients can be adjusted according to the actual network structure characteristics to balance the dominant role of core nodes and the local supporting role of neighboring nodes. By comparing the module affinity of conflicting nodes with each candidate module, they are reassigned to the module with the highest module affinity, thereby optimizing the assignment relationship of module boundaries and making the final SNP functional module partitioning result more accurate and reasonable.

[0125] Conflicting nodes are reassigned to the SNP functional modules with the highest module affinity, while simultaneously being removed from other SNP functional modules. After the reassignment of all conflicting nodes, each node ultimately belongs to only one SNP functional module, resulting in clearer and more reasonable boundaries between SNP functional modules, thus yielding the final SNP functional module partitioning result.

[0126] S43. For each SNP functional module, the association significance between the joint effect of all SNPs and the corrected unbiased phenotypic value is calculated using the multi-label association analysis method to obtain the association value of each SNP functional module, and the SNP functional modules that are significantly associated are selected according to the preset significance threshold.

[0127] In a preferred embodiment, the step of calculating the significance of the association between the joint effect of all SNPs and the corrected unbiased phenotypic value using a multi-label association analysis method for each SNP functional module, obtaining the association value of each SNP functional module, and selecting the SNP functional modules with significant association according to a preset significance threshold includes the following steps:

[0128] S431. Extract the genotype matrix of all SNP marker data in each SNP functional module and match it with the corrected unbiased phenotypic values ​​to construct the association analysis dataset for each SNP functional module.

[0129] In a preferred embodiment, the step of extracting the genotype matrix of all SNP marker data within each SNP functional module and matching it with the corrected unbiased phenotypic values ​​to construct the association analysis dataset for each SNP functional module includes the following steps:

[0130] S4311. For each SNP functional module, extract the genotype information of all SNP marker data in each individual of the breeding population, and convert the genotypes into numerical data using numerical encoding to construct an initial genotype matrix.

[0131] It should be noted that before performing association analysis on SNP functional modules, the genotype information of all SNPs within a module needs to be converted into a numerical data format suitable for statistical analysis. Specifically, for each SNP functional module, all SNP markers contained within that module are traversed, and for each individual in the breeding population, the genotype information at these SNP loci is extracted. Genotype information is typically represented in the form of allele combinations, such as AA, AB, and BB.

[0132] To perform subsequent mathematical calculations and statistical analysis, genotypes need to be quantified using numerical coding. Common coding methods include: for dimorphic SNPs, encoding the genotype as 0, 1, and 2 to represent the homozygous reference allele, heterozygous, and homozygous variant allele, respectively; or using dominant coding as 0 and 1; or using allelic dosage coding. The coding values ​​of each SNP in each individual are arranged in rows (individuals) and columns (SNPs) to form a two-dimensional matrix, i.e., the initial genotype matrix. The rows of this matrix correspond to the individuals in the breeding population, the columns correspond to the SNP markers within the modules, and the matrix elements are the numerical genotype codes of the corresponding individuals at the corresponding SNP loci.

[0133] S4312. Perform quality control on the initial genotype matrix by removing SNP markers with genotype deletion rates exceeding a preset threshold and filling the remaining missing values ​​with the mean to obtain a complete genotype matrix.

[0134] It should be noted that the initial genotype matrix may contain missing genotype values ​​due to insufficient sequencing depth, poor alignment quality, or some SNP markers with excessively high deletion rates in the population. These factors can interfere with subsequent association analyses, thus requiring quality control. Specifically, the genotype deletion rate for each SNP marker in the population is calculated, which is the proportion of individuals with a missing genotype at that SNP locus to the total number of individuals. A preset deletion rate threshold is set, typically 0.1 or 0.2 (i.e., 10% or 20%). SNP markers with deletion rates exceeding this threshold are removed from the matrix. For SNP markers that remain but have sporadic missing values, missing value imputation is performed. After removing high-deletion-rate SNPs and imputing remaining missing values, a complete genotype matrix is ​​obtained. This matrix no longer contains missing values, and all elements are numerical data.

[0135] Common imputation methods include mean imputation, which replaces missing values ​​with the arithmetic mean of all known genotypes encoded by the SNP locus. Mean imputation maintains the overall mean level of the SNP, is simple to calculate, and has minimal impact on subsequent analyses.

[0136] S4313. Match the complete genotype matrix with the corrected unbiased phenotypic values ​​to form the association analysis dataset for each SNP functional module.

[0137] It should be noted that the processed complete genotype matrix contains the genotype information of each individual at each SNP within each module, while the obtained corrected unbiased phenotypic values ​​contain the target trait phenotypic information for each individual. For association analysis, these two sets of data need to be matched one-to-one on an individual basis.

[0138] Specifically, ensure that the rows (individuals) in the genotype matrix are in exactly the same order as the individuals in the unbiased phenotypic values, or align them using individual IDs. For each SNP functional module, merge its corresponding complete genotype matrix with the unbiased phenotypic values ​​of the individuals corresponding to that module to form a complete dataset. This dataset contains the genotype information of all SNPs within that module (as independent variables), the phenotypic values ​​of each individual (as dependent variables), and possible covariate information. This dataset is the association analysis dataset for that SNP functional module.

[0139] S432. Calculate the significance of the association between the joint effect of all SNPs within each SNP functional module and the unbiased phenotypic value using the multi-label association analysis method, and obtain the association value of each SNP functional module.

[0140] In a preferred embodiment, the step of calculating the significance of the association between the joint effect of all SNPs within each SNP functional module and the unbiased phenotypic value using multi-label association analysis to obtain the association value for each SNP functional module includes the following steps:

[0141] S4321. Based on the association analysis dataset of each SNP functional module, select a multi-marker statistical model suitable for association analysis of high-dimensional genetic data from the preset multi-marker statistical model library.

[0142] It should be noted that when conducting association analysis between SNP functional modules and phenotypic traits, traditional single-marker association analysis methods cannot effectively capture the joint effects of SNPs within a module, as each SNP functional module contains multiple SNP markers and there may be complex linkage disequilibrium relationships and epistatic effects between SNPs. This can easily lead to the omission of signals from the combined effects of multiple minor SNPs. Therefore, it is necessary to select a statistical model specifically designed for multi-marker joint analysis. Based on the specific characteristics of the association analysis dataset for each SNP functional module, such as the number of individuals in the breeding population, the number of SNPs within the SNP functional module, the type of phenotypic data, and the linkage disequilibrium structure between SNPs, the most suitable model should be selected from a pre-defined multi-marker statistical model library.

[0143] The multi-label statistical model library includes sequence kernel association test model, multiple linear regression model, and Bayesian sparse linear mixture model.

[0144] The sequence kernel association test model employs a statistical structure based on variance component testing. By constructing a kernel function, it maps the genotypic information of multiple SNPs into a high-dimensional space, measuring the genetic similarity between individuals within a module. This similarity matrix is ​​then incorporated into the mixed linear model framework as the covariance structure of random effects. During model training, given unbiased phenotypic values ​​and covariates, the variance components are estimated using the restricted maximum likelihood method. A score test is then used to construct a test statistic, which follows a mixed chi-square distribution. The p-value of the overall module association is calculated using numerical integration or moment matching methods, thus achieving the goal of detecting the joint effect of modules without estimating the individual effects of each SNP.

[0145] The multiple linear regression model employs a linear combination regression structure, using the numerical genotype codes of all SNPs within the module as multiple independent variables to establish a linear relationship with the unbiased phenotypic values. The model form is y = Xβ + ε, where y is the phenotypic vector, X is a matrix containing multiple SNP genotypes, β is the effect coefficient vector of each SNP, and ε is the residual term. During model training, the least squares method or the maximum likelihood method is used for parameter estimation. The effect coefficients and standard errors of each SNP are obtained by solving the normal equation, and the significance of the joint effect of all SNPs is calculated based on the F-test or the likelihood ratio test. This model has a transparent structure, interpretable parameters, and the contribution direction and magnitude of each SNP can be directly obtained.

[0146] The Bayesian sparse linear mixture model employs a hierarchical Bayesian structure with sparse priors. It assumes that only a few SNPs have non-zero effects and shrinks the effect coefficients by introducing prior distributions. The Bayesian sparse linear mixture model takes the form y = Xβ + Zu + ε, where β is the sparse SNP effect vector and u is the polygenic background effect. During model training, Markov chain Monte Carlo (MCMC) sampling or variational Bayesian inference is used for posterior estimation. The posterior distribution of the effect coefficients of each SNP is obtained through iterative sampling. The posterior probability of each SNP being included in the model is calculated, and the overall association significance of the module is evaluated based on the posterior inclusion probability or Bayesian factors.

[0147] S4322. Based on the selected multi-label statistical model, calculate the joint effect statistic between all SNP labels and unbiased phenotypic values ​​in the SNP functional module and their corresponding original association values.

[0148] It should be noted that after selecting a multi-label statistical model, the association analysis dataset of the SNP functional module is used as input to run the selected model for parameter estimation and statistical inference. Different models will produce different joint effect statistics and association value representations: for the sequence kernel association test model, the statistic is the variance component test statistic, which approximately follows a mixed chi-square distribution, and the corresponding raw p-value, i.e., the raw association value, is obtained by calculation; for the multiple linear regression model, the statistic is the F-statistic or likelihood ratio statistic, and the raw p-value is calculated through the F-distribution or chi-square distribution; for the Bayesian model, the posterior probability of each SNP is obtained through Markov chain Monte Carlo sampling, and then the Bayesian factor or posterior inclusion probability of the module as a whole is calculated as a measure of association strength. Regardless of the model used, a raw association value is ultimately obtained to measure the significance of the association between all SNPs in the module as a whole and the phenotype. This raw association value reflects the degree of joint contribution of the SNP functional module to the target trait in terms of genetic structure and is the basic statistic for judging whether the SNP functional module significantly affects the phenotype.

[0149] S4323. Perform multiple hypothesis testing to correct the original association values, and use the corrected original association values ​​as the final association values ​​for each SNP functional module.

[0150] It should be noted that since multiple SNP functional modules are obtained, each SNP functional module needs to undergo the aforementioned association analysis. This results in a set of raw association values ​​corresponding to multiple SNP functional modules, such as multiple p-values. Performing so many statistical tests simultaneously significantly increases the risk of false positives, meaning that modules that are actually unrelated to the phenotype may be incorrectly classified as significantly associated. To control the accumulation of Type I errors caused by this multiple testing, multiple hypothesis testing correction is needed for the raw association values ​​of all SNP functional modules. Commonly used correction methods include false positive rate control (FPR) methods. This method balances false positives and statistical power by controlling the proportion of false positives, and is suitable for large-scale screening scenarios where a certain proportion of false positives is tolerated while detecting true signals. It is widely used in genetic association analysis.

[0151] S433. Compare the correlation values ​​of SNP functional modules with a preset significance threshold, and select SNP functional modules with correlation values ​​less than the significance threshold as significantly correlated SNP functional modules.

[0152] Specifically, after obtaining the final association value of each SNP functional module after multiple hypothesis testing correction, these modules need to be screened according to a pre-set significance threshold to determine which modules are statistically significantly associated with the target trait. The pre-set significance threshold is usually determined based on the research purpose and the rigor of the multiple tests. For example, in genome association studies, FDR < 0.05 is often used as the significance threshold, i.e., controlling the false discovery rate to within 5%.

[0153] S44. For each significantly associated SNP functional module, and using the transformation of the cumulative distribution function, select several SNPs as representative markers in each SNP functional module, and construct the screened whole genome SNP marker data from all representative markers.

[0154] It should be noted that after obtaining the significantly associated SNP functional modules, each SNP functional module may contain dozens or even hundreds of SNP markers, and these SNPs have varying degrees of functional association within the module. In order to further reduce marker redundancy and extract the most representative SNPs for subsequent breeding practices, it is necessary to rank and screen the SNPs within each significantly associated module according to their importance.

[0155] Specifically, for each significantly associated SNP functional module, a module-wide network is constructed, consisting of sub-networks with all SNPs within the module as nodes and their biofunctional correlation coefficients as edge weights. Simultaneously, an individual network is constructed for each SNP within the module, consisting only of the SNP itself and its connections to other SNPs within the module. The cumulative distribution function is used to transform the edge weight distributions in both the module-wide network and each SNP's individual network into corresponding cumulative distribution curves. Then, the Kolmogorov-Smirnov statistic (KS distance) is used to calculate the distance between the cumulative distribution of each SNP's individual network and the cumulative distribution of the module-wide network. A smaller KS distance indicates a more similar local network structure to the overall module structure, meaning the SNP is more representative of the module and better reflects its functional characteristics. All SNPs within a module are sorted in ascending order of KS distance. A few top-ranked SNPs from each module are selected, for example, the top 5, top 10, or a certain percentage (e.g., the top 20%), as representative markers for that SNP functional module. These representative markers from all significantly associated SNP functional modules are then merged to form the final, filtered genome-wide SNP marker dataset.

[0156] S5. Using gene set and transcript information, functional annotation was performed on the screened whole-genome SNP marker data to obtain molecular markers for the breeding of Snowflake Black Cattle.

[0157] It is important to note that each SNP in the selected whole-genome SNP marker dataset requires detailed functional annotation. Specifically, using the obtained gene set and transcript information, combined with public bioinformatics databases, each selected SNP marker is subjected to multi-dimensional functional annotation: this includes determining the physical location of the SNP on the reference genome, identifying whether it is located within a gene, between genes, or in a regulatory region; if the SNP is located within a gene, using transcript information to determine whether it causes amino acid changes or affects splicing; then, combining eQTL analysis results, annotating whether the SNP regulates the expression of other genes and its target genes and the direction of effect; based on gene set information, annotating the biological pathways, molecular functions, and cellular components involved by the genes in which the SNP is located or regulates; and searching existing literature and databases to understand the progress of research on the association between the SNP or its gene and meat quality traits in other species or varieties. Through the above multi-dimensional functional annotation, a complete functional description of each SNP is finally formed, including its location information, variation type, regulatory relationship, functional pathway, and relevant literature support. Functionally annotated SNP markers are molecular markers that can be used for the breeding of Snowflake Black Cattle. They can be directly applied to subsequent marker-assisted selection or genomic selection breeding practices, providing precise molecular tools for the genetic improvement of Snowflake Black Cattle.

[0158] like Figure 2 As shown, according to another embodiment of the present invention, a Snowflake Black Cattle breeding system based on genome-wide association analysis is provided, the system comprising:

[0159] Data acquisition module 1 is used to collect blood and tissue samples from the Snowflake Black Cattle breeding population, and extract PacBio long-read HiFi sequencing data and second-generation whole genome resequencing data from the blood samples, as well as extract Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data from the tissue samples.

[0160] Assembly sequencing module 2 is used to assemble the genome based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data to obtain a chromosome-level reference genome sequence; and to annotate the reference genome sequence with gene structure based on multi-tissue full-length transcriptome sequencing data to obtain gene set and transcript information.

[0161] Data alignment module 3 is used to align second-generation whole-genome resequencing data to the reference genome sequence and detect the whole-genome SNP marker dataset of the breeding population;

[0162] Data filtering module 4 is used to collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and combined with the whole genome SNP marker dataset, to filter whole genome SNP marker data that meet the preset significance threshold through whole genome association analysis;

[0163] Functional annotation module 5 is used to perform functional annotation on the screened whole-genome SNP marker data using gene set and transcript information to obtain molecular markers for the breeding of Snowflake Black Cattle.

[0164] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, optical storage, etc.) containing computer-usable program code.

[0165] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A breeding method for Snowflake Black Cattle based on genome-wide association analysis, characterized in that, The method includes the following steps: S1. Collect blood and tissue samples from the Snowflake Black Cattle breeding population, and extract PacBio long-read HiFi sequencing data and second-generation whole-genome resequencing data from the blood samples, as well as Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data from the tissue samples. S2. Based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data, genome assembly was performed to obtain a chromosome-level reference genome sequence; based on multi-tissue full-length transcriptome sequencing data, gene structure annotation was performed on the reference genome sequence to obtain gene set and transcript information; S3. Align the second-generation whole-genome resequencing data to the reference genome sequence and detect the whole-genome SNP marker dataset of the breeding population; S4. Collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and combine them with the whole genome SNP marker dataset. Through whole genome association analysis, screen whole genome SNP marker data that meet the preset significance threshold. S5. Using gene set and transcript information, functional annotation was performed on the screened whole-genome SNP marker data to obtain molecular markers for the breeding of Snowflake Black Cattle.

2. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 1, characterized in that, The process of collecting standardized phenotypic data from all individuals in the Snowflake Black Cattle breeding population, and combining this data with a genome-wide SNP marker dataset, to screen genome-wide SNP marker data that meet a preset significance threshold through genome-wide association analysis includes the following steps: S41. Collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and obtain corrected unbiased phenotypic values ​​by correcting the standardized phenotypic data. S42. Based on the whole genome SNP marker dataset, construct an SNP relationship network, and use local co-optimization to divide the SNP relationship network to obtain several SNP functional modules. S43. For each SNP functional module, the association significance between the joint effect of all SNPs and the corrected unbiased phenotypic value is calculated using the multi-label association analysis method to obtain the association value of each SNP functional module, and the SNP functional modules that are significantly associated are selected according to the preset significance threshold. S44. For each significantly associated SNP functional module, and using the transformation of the cumulative distribution function, select several SNPs as representative markers in each SNP functional module, and construct the screened whole genome SNP marker data from all representative markers.

3. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 2, characterized in that, The process of constructing an SNP relationship network based on a whole-genome SNP marker dataset and dividing the SNP relationship network into several SNP functional modules using a local co-optimization approach includes the following steps: S421. Based on the whole genome SNP marker dataset, calculate the biological functional association coefficient between any two SNP marker datasets; S422. Construct an SNP relationship network with SNP-labeled data as nodes and biological function correlation coefficients as edge weights; S423. By calculating the validity of each node in the SNP relationship network and combining it with local collaborative optimization, the SNP relationship network is divided into several SNP functional modules.

4. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 3, characterized in that, The calculation of the biological functional association coefficient between any two SNP marker data based on the whole genome SNP marker dataset includes the following steps: S4211. Based on the Hi-C chromosome conformation, a chromosome-level chromatin interaction matrix is ​​constructed using the captured sequencing data, and the normalized contact frequency between any two SNP marker data regions is extracted as the spatial interaction strength coefficient. S4212. Based on gene set and transcript information, combined with a pre-set biological pathway database, pathway annotation is performed on each SNP marker data, and semantic similarity is used to calculate the functional similarity coefficient between any two SNP marker data. S4213. Perform quantitative trait locus localization analysis on multi-tissue full-length transcriptome sequencing data to obtain the target genes associated with each SNP as a regulatory site and their effect direction and intensity, and calculate the correlation coefficient of any two SNP marker data on the regulatory expression profile based on the effect direction and intensity. S4214. The spatial interaction strength coefficient, functional similarity coefficient, and correlation coefficient are weighted and summed according to preset weight coefficients to obtain the biological functional correlation coefficient.

5. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 3, characterized in that, The process of dividing the SNP relationship network into several SNP functional modules by calculating the effectiveness of each node in the SNP relationship network and combining it with local collaborative optimization includes the following steps: S4231. Mark all nodes as unvisited and initialize the functional module list to empty; S4232. Calculate the validity of each node in the SNP relationship network based on the sum of the biological functional correlation coefficients of all edges connected to the node, and select the node with the highest validity as the core node. S4233. Starting from the core node, perform local collaborative expansion processing to generate candidate SNP functional modules, and mark the core node as visited. S4234. Repeat steps S4232 to S4233 until all nodes are marked as visited, and obtain the preliminary SNP functional module partitioning results. S4235. Analyze the conflict nodes in the preliminary SNP functional module partitioning results, and based on the analysis results, reallocate the preliminary SNP functional module partitioning results to obtain the final SNP functional module partitioning results.

6. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 5, characterized in that, The process of performing local collaborative expansion starting from the core node to generate candidate SNP functional modules and marking the core node as visited includes the following steps: S42331. Based on the biological functional correlation coefficient between the core node and its neighboring nodes, calculate the correlation strength between the core node and all its neighboring nodes, and sort the neighboring nodes in descending order of correlation strength. S42332. Initialize the candidate SNP function module, mark the core node as visited and include it in the candidate module, select the first-ranked neighbor node to add to the candidate SNP function module and mark it as visited, integrate the neighbor nodes of all nodes in the candidate SNP function module to form a set of neighbors to be expanded, and remove visited nodes at the same time. S42333. Calculate the internal and external degrees of the current candidate SNP functional modules, and calculate the fitness of each node in the neighbor set to be expanded based on the internal and external degrees. S42334. Add neighbor nodes with fitness greater than zero to the candidate SNP function module and mark them as visited. If no node in the neighbor set to be expanded satisfies fitness greater than zero, end the local collaborative expansion process.

7. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 2, characterized in that, For each SNP functional module, the association significance between the joint effect of all SNPs and the corrected unbiased phenotypic value is calculated using multi-label association analysis to obtain the association value of each SNP functional module. The steps for selecting significantly associated SNP functional modules based on a preset significance threshold are as follows: S431. Extract the genotype matrix of all SNP marker data in each SNP functional module and match it with the corrected unbiased phenotypic values ​​to construct the association analysis dataset for each SNP functional module. S432. Calculate the significance of the association between the joint effect of all SNPs within each SNP functional module and the unbiased phenotypic value using the multi-label association analysis method, and obtain the association value of each SNP functional module. S433. Compare the correlation values ​​of SNP functional modules with a preset significance threshold, and select SNP functional modules with correlation values ​​less than the significance threshold as significantly correlated SNP functional modules.

8. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 7, characterized in that, The steps involved in extracting the genotype matrix of all SNP marker data within each SNP functional module and matching it with the corrected unbiased phenotypic values ​​to construct the association analysis dataset for each SNP functional module are as follows: S4311. For each SNP functional module, extract the genotype information of all SNP marker data in each individual of the breeding population, and convert the genotypes into numerical data using numerical encoding to construct an initial genotype matrix. S4312. Perform quality control on the initial genotype matrix by removing SNP markers with genotype deletion rates exceeding a preset threshold and filling the remaining missing values ​​with the mean to obtain a complete genotype matrix. S4313. Match the complete genotype matrix with the corrected unbiased phenotypic values ​​to form the association analysis dataset for each SNP functional module.

9. The method for breeding Snowflake Black Cattle based on genome-wide association analysis according to claim 7, characterized in that, The step of calculating the significance of the association between the joint effect of all SNPs within each SNP functional module and the unbiased phenotypic value using multi-label association analysis to obtain the association value for each SNP functional module includes the following steps: S4321. Based on the association analysis dataset of each SNP functional module, select a multi-marker statistical model suitable for association analysis of high-dimensional genetic data from the preset multi-marker statistical model library. S4322. Based on the selected multi-label statistical model, calculate the joint effect statistic between all SNP labels and unbiased phenotypic values ​​in the SNP functional module and their corresponding original association values. S4323. Perform multiple hypothesis testing to correct the original association values, and use the corrected original association values ​​as the final association values ​​for each SNP functional module.

10. A Snowflake Black Cattle breeding system based on genome-wide association analysis, used to implement the Snowflake Black Cattle breeding method based on genome-wide association analysis as described in any one of claims 1-9, characterized in that, The system includes: The data acquisition module is used to collect blood and tissue samples from the Snowflake Black Cattle breeding population, and extract PacBio long-read HiFi sequencing data and second-generation whole-genome resequencing data from the blood samples, as well as extract Hi-C chromosome conformation capture sequencing data and multi-tissue full-length transcriptome sequencing data from the tissue samples. The assembly sequencing module is used to assemble the genome based on PacBio long-read HiFi sequencing data and Hi-C chromosome conformation capture sequencing data to obtain a chromosome-level reference genome sequence; and to annotate the reference genome sequence with gene structure based on multi-tissue full-length transcriptome sequencing data to obtain gene set and transcript information. The data alignment module is used to align second-generation whole-genome resequencing data to a reference genome sequence and to detect whole-genome SNP marker datasets of the breeding population. The data filtering module is used to collect standardized phenotypic data of all individuals in the Snowflake Black Cattle breeding population, and combined with the whole genome SNP marker dataset, to filter whole genome SNP marker data that meet the preset significance threshold through genome-wide association analysis; The functional annotation module is used to perform functional annotation on the screened whole-genome SNP marker data using gene set and transcript information, so as to obtain molecular markers for the breeding of Snowflake Black Cattle.

Images