Guiding method and system for detecting animal inbreeding level based on whole genome sequencing method

Abstract

The invention provides a guiding method and system for detecting animal inbreeding level based on a whole genome sequencing method, and belongs to the field of genomics and molecular biology. The problem that errors are large when ROH detection is carried out on low-quality data is solved. The method comprises the following steps: acquiring real reference genome data and whole genome sequencing data, and performing preprocessing, namely screening autosomes without base deletion sequences to construct a new reference genome, calculating real inbreeding level parameters and generating a haplotype genome file; the simulation whole genome sequencing data set comprises a sequencing depth gradient distribution data set, a read length gradient distribution data set and a reference genome continuity gradient distribution data set; respectively carrying out homozygous fragment detection by using the simulated whole genome sequencing data set so as to evaluate the influence of data quality parameters on the detection accuracy; and optimizing homozygous fragment detection software parameters based on an evaluation result. The method is used in the field of inbreeding evaluation of endangered species and small populations.

Description

Technical Field

[0001] This invention belongs to the field of genomics and molecular biology, and in particular relates to a guiding method for detecting the level of inbreeding in animals based on whole-genome sequencing. Background Technology

[0002] Inbreeding is often accompanied by inbreeding depression, posing a serious threat to the long-term survival of small, isolated populations, especially endangered species. The genetic characterization of genome-wide inbreeding is a key indicator for assessing the extinction risk of endangered species and small populations. When inbreeding occurs, homologous chromosome segments (IBDs) are passed on to offspring, leading to homozygosity in the offspring genome. Currently, these IBD segments can be accurately characterized by identifying homozygous regions throughout the genome using whole-genome sequencing (WGS) data. Homozygous segments (ROHs) in the whole genome serve as a core indicator for assessing the level of genome-wide inbreeding, accurately estimating the past and present inbreeding status of a population without relying on pedigree records. It is generally believed that the ROH segments from the first generation of inbreeding are the longest. After several generations of continuous inbreeding, ROH segments may be broken by meiotic recombination. Therefore, longer ROH segments (>1 Mb in mammals) typically reflect recent inbreeding events (within 50 generations), while shorter ROH segments reflect historical inbreeding events. In addition, the ratio F of the total length of the ROH segment on autosomes to the total physical length of autosomes is usually used. ROH To quantify the degree of inbreeding within a population, F ROH The higher the value, the more severe the inbreeding in the population. Currently, the length distribution and proportion of ROH segments in the genome are widely used to assess recent and historical inbreeding events in a population, and further to assess the extinction risk of endangered species or prevent population decline during the breeding of economically important animals.

[0003] ROH fragments can typically reach millions of base pairs in length, making accurate ROH detection highly dependent on high-quality reference genomes and resequencing data. In practice, sampling endangered species is extremely difficult; their DNA samples are often obtained through non-endogenous methods such as feces and hair, resulting in severely fragmented DNA and low sequencing depth. Such low-quality data can introduce significant errors in subsequent ROH detection. Furthermore, even using the same sequencing data, different contiguous reference genomes will result in varying ROH fragment lengths and F... ROHInconsistencies can also occur. Furthermore, some endangered species present challenges in obtaining high-quality DNA for sampling, making it impossible to assemble a reference genome. In such cases, it's necessary to use a reference genome from a closely related species for ROH detection. If the target population is genetically distant from the reference genome, inaccurate ROH fragment detection may occur, especially for long ROH fragments resulting from recent inbreeding, which may be incorrectly broken into many smaller ROH fragments. PLINK software is one of the most commonly used tools for ROH detection based on WGS data; however, when data quality is low, its default parameters may not be suitable, easily producing false positive or false negative ROH fragments, rendering inbreeding assessment results unreliable.

[0004] Currently, there are no standard guidelines for the data quality requirements of ROH testing. Therefore, it is urgent to establish a set of standard guidelines based on whole-genome sequencing technology to accurately detect the inbreeding level of animals, clarify the key thresholds of the reference genome and WGS data quality parameters used for ROH testing, and optimize various parameters of the ROH testing software for low-quality data, thereby solving the above problems and improving the accuracy of ROH testing. Summary of the Invention

[0005] Therefore, this invention aims to propose a standard guidance method for detecting the level of inbreeding in animals based on whole-genome sequencing, in order to solve the problem of large errors when using low-quality data for ROH detection.

[0006] To achieve the above objectives, the present invention employs the following technical solution: a guided method for detecting animal inbreeding levels based on whole-genome sequencing, the method comprising: Step S1: Obtain real high-quality reference genome data and high-quality whole genome sequencing data of a real population from a public database, and preprocess the obtained data. The preprocessing includes screening complete autosomal sequences with no base deletions from telomeres to telomeres from the high-quality reference genome to construct a new reference genome, and using the new genome as a reference, performing variant detection based on high-quality resequencing data to identify the genetic variations in the genome of each individual in the real population. Step S2: Calculate the inbreeding evaluation parameters ROH and F for different individuals in the real population. ROH The variation data of each individual will be restored to the reference genome, and two haplotype genome files with real variation characteristics will be generated for each individual. Step S3: Simulate the reference genome dataset and the whole genome sequencing dataset using the preprocessed real data. The simulated whole genome sequencing dataset includes a dataset with sequencing depth gradient distribution, a dataset with sequencing read length gradient distribution, and a dataset with reference genome continuity gradient distribution. Step S4: Homozygous fragment detection was performed using simulated reference genome datasets and whole genome sequencing datasets to assess the impact of data quality parameters on detection accuracy; Step S5: Optimize the parameters of the ROH fragment detection software based on the evaluation results for low-quality sequencing data.

[0007] Furthermore, the F ROH The calculation method is as follows:

[0008] Among them, F ROH This represents the ratio of the total physical length of autosomes. This represents the total length of the ROH segment on an autosome. The total length of the autosome is denoted as 100 kb; the minimum ROH fragment length threshold is set to 100 kb.

[0009] Furthermore, a preferred method is proposed, wherein the sequencing depth gradient includes: 100×, 60×, 30×, 20×, 15×, 12×, 10×, 8×, 6×, 5×, 4×, 3.5×, 3×, 2.5×, 2×, and 1×; the sequencing read length gradient includes: 30bp, 50bp, 75bp, 100bp, and 150bp; the continuity of the reference genome is characterized by the value of contig N50, and the gradient in the reference genome includes: 25.91Kb, 44.90Kb, 88.84Kb, 177.75Kb, 304.10Kb, 600.66Kb, 1.20Mb, 2.05Mb, 4.06Mb, 8.13Mb, 13.89Mb, 27.38Mb, 54.88Mb, and 93.60Mb.

[0010] Furthermore, a preferred method is proposed: in step S3, the data in the simulated whole-genome sequencing dataset is generated using wgsim software with parameters set to fixed read length and without introducing sequencing errors; and the haplotype genome file is generated by replacing variant sites in the reference genome. When the variant site is a 1 / 1 genotype, the bases at that position for both haplotypes are replaced with the allele, while for the 0 / 1 genotype, only one haplotype is replaced.

[0011] Furthermore, a preferred method is proposed, wherein the homozygous fragment detection in step S4 includes: Using BWA software mem The algorithm compares whole-genome sequencing data with reference gene sequencing data; The generated BAM files are preprocessed, including sorting and deduplication. Population-level variation detection was performed on the preprocessed BAM files, including; The HaplotypeCaller algorithm from the GATK toolkit is used to call the original mutations of each individual and generate a gVCF file; Use the CombineGVCFs tool in the GATK toolkit to merge these gVCF files into a single gVCF file, and use the GenotypeGVCFs toolkit in the GATK toolkit to perform joint mutation detection and generate a population-level VCF file. Single nucleotide polymorphisms (SNPs) were extracted from all variants and hard-filtered. The SNPs were further filtered using VCFtools software to remove all non-biallelic SNPs and SNPs with a deletion rate of more than 20%. ROH detection was performed using PLINK software, and the VCF file was converted to PLINK file format; the --indep-pairwise command was used for linkage imbalance LD filtering, with parameters set to a window size of 50Kb and a correlation index. The threshold is 0.9; Use the `--homozyg` command to perform homozygous fragment detection, including: Each sliding window contains at least 20 SNPs, with the SNP density within the window being one out of every 50 SNPs, and the physical length of the ROH fragment is at least 100 kb.

[0012] Furthermore, a preferred embodiment is proposed, wherein step S4 further includes using F ROH Error rate assessment at different sequencing depths F ROH The accuracy of the detection, F ROH The formula for calculating the error rate is: .

[0013] Furthermore, a preferred embodiment is proposed, wherein step S5 includes: The command `plink --file filename --indep-pairwise n1 n2 n3` is used to calculate the linkage disequilibrium LD between SNP pairs within a window of n1, and the process is performed by advancing n2 SNPs in each iteration; if the correlation index of any SNP pair within the window is... If the value exceeds n3, mark one of the SNPs as redundant and delete it; continuously adjust the window size n1 of the command --indep-pairwise from 1Mb to 5Kb. The threshold n3 was continuously adjusted from 0.01 to 0.99; ROH detection is performed using the commands `--homozyg --homozyg-window-snp n4 --homozyg-density n5 --homozyg-kb n6 --homozyg-snp n7`. This step scans the genome using a sliding window of n4 SNPs, allowing a maximum of n5 SNPs per Mb region to control SNP density. Regions longer than n6 and containing at least n7 SNPs are identified as ROHs. The parameter n4 of the command `--homozyg-window-snp` is continuously adjusted from 100 SNPs to 5 SNPs, and the parameter n7 of the command `--homozyg-snp` is continuously adjusted from 50 SNPs to 300 SNPs.

[0014] Based on the same inventive concept, this invention also proposes a guidance system for detecting animal inbreeding levels based on whole-genome sequencing, the system comprising: The preprocessing unit is used to obtain real high-quality reference genome data and high-quality whole genome sequencing data of a real population from a public database, and to preprocess the obtained data. The preprocessing includes screening complete autosomal sequences with no base deletions from telomeres to telomeres from the high-quality reference genome to construct a new reference genome, and using the new genome as a reference, performing variant detection based on high-quality resequencing data to identify the genetic variations in the genome of each individual in the real population. The haplotype genome file generation unit is used to calculate the inbreeding evaluation parameters ROH and F for different individuals in a real population. ROH The variation data of each individual will be restored to the reference genome, and two haplotype genome files with real variation characteristics will be generated for each individual. The simulation unit is used to simulate a reference genome dataset and a whole genome sequencing dataset using preprocessed real data. The simulated whole genome sequencing dataset includes a dataset with sequencing depth gradient distribution, a dataset with sequencing read length gradient distribution, and a dataset with reference genome continuous gradient distribution. The homozygous fragment detection unit is used to perform homozygous fragment detection using simulated reference genome datasets and whole genome sequencing datasets, respectively, to assess the impact of data quality parameters on detection accuracy. The optimization unit is used to optimize the parameters of the ROH fragment detection software based on the evaluation results.

[0015] Based on the same inventive concept, the present invention also proposes a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the processor runs the computer program stored in the memory, the processor executes a guided method for detecting animal inbreeding levels based on whole-genome sequencing according to any one of the above.

[0016] Based on the same inventive concept, the present invention also proposes a computer-readable storage medium storing a computer program that, when executed by a processor, performs the steps of the guided method for detecting animal inbreeding levels based on whole-genome sequencing as described above.

[0017] Compared with the prior art, the beneficial effects of the present invention are: Existing technologies for ROH detection lack clear standards regarding the required quality of reference genome and sequencing data, leading to unpredictable results. This invention pioneers a method for systematic testing using simulated gradient datasets. By controlling variables, it simulates reference genome data with varying continuity, whole-genome sequencing data with different sequencing depths and read lengths, and accurately assesses the accuracy of ROH detection. By providing clear data quality thresholds and optimized parameter schemes, this invention effectively reduces F (Failure Rate) errors caused by poor data quality. ROH Overestimation or underestimation. Experimental data in the disclosure document indicate that, at the recommended sequencing depth of 15x, F... ROH The error rate can be reduced to an extremely low level, namely 0.84% ​​± 1.71%, ensuring the accuracy of inbreeding level assessment.

[0018] Existing PLINK technology uses a one-size-fits-all approach with default parameter settings, which cannot adapt to fluctuations in data quality. This invention discovers an intrinsic correlation between sequencing depth and reference genome continuity, and establishes a dynamic parameter optimization mechanism based on this. For example, it was found that for low-depth data, adjusting the combination of parameters such as "--indep-pairwise" and "--homozyg-window-snp" significantly corrects the estimation bias of FROH.

[0019] To address different needs such as high-quality samples, low-quality samples (e.g., fecal samples), and highly inbred individuals, this invention provides a simple and easy-to-use parameter scheme. In particular, it allows for parameter adjustments for situations where data quality may be poor due to the difficulty in sampling endangered species. This helps researchers design sequencing schemes in advance based on actual conditions, avoid repeated sequencing, and reduce experimental and time costs.

[0020] This invention is applicable to the accurate assessment of inbreeding in the conservation of endangered species and small populations, as well as in the breeding of economically important animals. Attached Figure Description

[0021] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 This is a flowchart illustrating the ROH detection process under different conditions as described in this invention; Figure 2 This is a schematic diagram illustrating the total number of ROH fragments detected using reference genomes with different continuity as described in this invention; Figure 3 F is calculated for the 8 sets of simulated WGS data described in this invention. ROH The trends at different sequencing depths are shown in the figure, where the numbers represent the true F. ROH value; Figure 4 The F-value for ROH detection calculated using simulated WGS data at different sequencing depths as described in this invention. ROH Error rate; Figure 5 This is a schematic diagram illustrating the differences in ROH detection results using simulated WGS data with different sequencing read lengths as described in this invention; Figure 6 This is a schematic diagram comparing the FROH results calculated using real WGS data of the same subspecies and reference genomes at different genetic distances, as described in this invention, and a schematic diagram comparing the number of ROH fragments larger than 5Mb detected. Figure 7 This is a flowchart illustrating the accurate detection of ROH based on whole-genome sequencing technology in protective genomics as described in this invention. Detailed Implementation

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. It should be noted that, unless otherwise specified, the embodiments and features in the embodiments of the present invention can be combined with each other, and the described embodiments are only some embodiments of the present invention, not all embodiments.

[0023] Implementation Method 1, see Figure 1 This embodiment describes a guided method for detecting inbreeding levels in animals based on whole-genome sequencing. The method includes: Step S1: Obtain real high-quality reference genome data and high-quality whole genome sequencing data of a real population from a public database, and preprocess the obtained data. The preprocessing includes screening complete autosomal sequences with no base deletions from telomeres to telomeres from the high-quality reference genome to construct a new reference genome, and using the new genome as a reference, performing variant detection based on high-quality resequencing data to identify the genetic variations in the genome of each individual in the real population. Step S2: Calculate the inbreeding evaluation parameters ROH and F for different individuals in the real population. ROH The variation data of each individual will be restored to the reference genome, and two haplotype genome files with real variation characteristics will be generated for each individual. Step S3: Simulate the reference genome dataset and the whole genome sequencing dataset using the preprocessed real data. The simulated whole genome sequencing dataset includes a dataset with sequencing depth gradient distribution, a dataset with sequencing read length gradient distribution, and a dataset with reference genome continuity gradient distribution. Step S4: Homozygous fragment detection was performed using simulated reference genome datasets and whole genome sequencing datasets to assess the impact of data quality parameters on detection accuracy; Step S4, homozygous fragment detection, includes: Using BWA software mem The algorithm compares whole-genome sequencing data with reference gene sequencing data; The generated BAM files are preprocessed, including sorting and deduplication. Population-level variation detection was performed on the preprocessed BAM files, including; The HaplotypeCaller algorithm from the GATK toolkit is used to call the original mutations of each individual and generate a gVCF file; Use the CombineGVCFs tool in the GATK toolkit to merge these gVCF files into a single gVCF file, and use the GenotypeGVCFs toolkit in the GATK toolkit to perform joint mutation detection and generate a population-level VCF file. Single nucleotide polymorphisms (SNPs) were extracted from all variants and hard-filtered. The SNPs were further filtered using VCFtools software to remove all non-biallelic SNPs and SNPs with a deletion rate of more than 20%. ROH detection was performed using PLINK software, and the VCF file was converted to PLINK file format; the --indep-pairwise command was used for linkage imbalance LD filtering, with parameters set to a window size of 50Kb and a correlation index. The threshold is 0.9; Use the `--homozyg` command to perform homozygous fragment detection, including: Each sliding window contains at least 20 SNPs, with the SNP density within the window being one out of every 50 SNPs, and the physical length of the ROH fragment is at least 100Kb. Step S4 also includes using F ROH Error rate assessment at different sequencing depths F ROH The accuracy of the detection, F ROH The formula for calculating the error rate is: ; Step S5: Optimize the parameters of the ROH fragment detection software based on the evaluation results, including: The command `plink --file filename --indep-pairwise n1 n2 n3` is used to calculate the linkage disequilibrium LD between SNP pairs within a window of n1, and the process is performed by advancing n2 SNPs in each iteration; if the correlation index of any SNP pair within the window is... If the value exceeds n3, mark one of the SNPs as redundant and delete it; continuously adjust the window size n1 of the command --indep-pairwise from 1Mb to 5Kb. The threshold n3 was continuously adjusted from 0.01 to 0.99; ROH detection is performed using the commands `--homozyg --homozyg-window-snp n4 --homozyg-density n5 --homozyg-kb n6 --homozyg-snp n7`. This step scans the genome using a sliding window of n4 SNPs, allowing a maximum of n5 SNPs per Mb region to control SNP density. Regions longer than n6 Kb and containing at least n7 SNPs are identified as ROHs. The parameter n4 of the command `--homozyg-window-snp` is continuously adjusted from 100 SNPs to 5 SNPs, and the parameter n7 of the command `--homozyg-snp` is continuously adjusted from 50 SNPs to 300 SNPs.

[0024] Specifically, the method proposed in this embodiment includes: Real reference genome data and whole-genome sequencing (WGS) data were obtained from public databases and preprocessed. The preprocessed real data were used to simulate whole-genome sequencing data with fixed read lengths and determined inbreeding levels, and this data was used to extract whole-genome sequencing data at different depths. Simulated whole-genome sequencing data with different read lengths at fixed sequencing depths were also performed. Finally, reference genome data with different continuity were simulated. This implementation method proposes obtaining reference genome data and whole genome sequencing data from public databases and performing preprocessing, specifically including: Reference genome: Select a high-quality reference genome and screen autosomes without telomere deletion features from telomere to telomere to construct a new reference genome; Calculate the true inbreeding level parameter F ROH : Calculate the true F-value of population whole-genome sequencing data downloaded from public databases ROH Value, retain F ROH Individual WGS data with values ​​exhibiting a gradient distribution are used for subsequent simulation analysis; The true inbreeding level parameter F ROH The calculation method is as follows:

[0025] in, This represents the total length of the ROH fragment on the euchromosomal WGS body. The total length of the autosome is 100 kb; the ROH fragment length threshold is set to 100 kb.

[0026] Generate haplotype genome files for each individual's real data: Align each individual's real WGS data to a new reference genome for variant detection to obtain a VCF file. Replace the original sites in the new reference genome with the corresponding variant sites in the VCF file to obtain two haplotype genome files for each individual (e.g., Hap1.fasta, Hap2.fasta). Specifically, for genotypes displayed as "1 / 1" in the VCF file, the original sites are replaced with variant sites in both Hap1.fasta and Hap2.fasta; for genotypes of "0 / 1", site replacement is performed only in Hap1.fasta or Hap2.fasta.

[0027] In this embodiment, preprocessed real data is used to simulate a fixed read length, F. ROH The WGS data with defined values ​​specifically includes: Use haplotype genome files generated based on real data for each individual as input.

[0028] WGS data for determining the ROH ratio was simulated using wgsim software with parameters "-1 100 -2 100 -e 0".

[0029] Determined F ROH The value is limited to the true F calculated using real whole-genome sequencing data. ROH value.

[0030] The read length is fixed at 100bp paired-end reads and does not introduce sequencing errors.

[0031] Reference genome data of the same subspecies, the same genus, and the same family, as well as whole genome sequencing data of the same subspecies, are obtained from public databases as test data for the accuracy of ROH detection under different genetic distances between the reference genome and the target population.

[0032] ROH detection was performed using simulated whole-genome sequencing data with different sequencing depth gradients, fragment length gradients, and reference genome continuity gradients to evaluate optimal and minimum data quality parameters. Simulated WGS data were randomly selected, and ROH detection was performed using real data at different genetic distances to assess the limitations of using closely related species as reference genomes for ROH detection. Finally, the optimization of parameters when using low-quality data for ROH detection was evaluated. The quality of ROH detection results was assessed using F... ROH Error rate is used to characterize the error rate.

[0033] Randomly sampled simulated WGS data, specifically including: WGS data from simulations at different depths were randomly extracted using the Seqkit software. The different depths are limited to 100×, 60×, 30×, 20×, 15×, 12×, 10×, 8×, 6×, 5×, 4×, 3.5×, 3×, 2.5×, 2× and 1×.

[0034] Using preprocessed real data to simulate WGS data of different read lengths with fixed sequencing depth, specifically including: Use haplotype genome files generated based on real data for each individual as input.

[0035] F was determined using wgsim software simulation. ROH The WGS data values ​​are set with gradients of 30bp to 150bp using parameters "-1 and -2", and parameter "-e" is always limited to 0.

[0036] The sequencing depth was fixed at 20×, and the WGS data of different read lengths were limited to 30bp, 50bp, 75bp, 100bp and 150bp.

[0037] This implementation simulates different reference genome continuity gradients, including: The reference genome is characterized by the value of contig N50. A deletion sequence (N) is randomly inserted into the new reference genome to simulate a fragmented genome; the inserted deletion sequence is limited to 300 bp "N". Evaluate the contig N50 and longest contig length for each simulated reference genome; Thirteen simulated reference genomes with contig N50 values ​​ranging from 25.91 Kb to 93.60 Mb were selected for subsequent analysis; the gradients in the reference genomes included: 25.91 Kb, 44.90 Kb, 88.84 Kb, 177.75 Kb, 304.10 Kb, 600.66 Kb, 1.20 Mb, 2.05 Mb, 4.06 Mb, 8.13 Mb, 13.89 Mb, 27.38 Mb, 54.88 Mb, and 93.60 Mb.

[0038] In this embodiment, ROH detection is performed using simulated data, specifically including: The simulated data used for ROH detection excludes CNV (copy number variation) regions; ROH detection was performed using datasets with different sequencing depths, different sequencing read lengths, and different contig N50 reference genome datasets, with a unique control variable during detection. ROH detection was performed using datasets of different sequencing depths, with the most complete new reference genome used and the whole genome sequencing data read length limited to 100 bp. The simulated whole-genome sequencing data used for ROH detection with different sequencing read length datasets had a depth of 20×, and the reference genome contig N50s used were 93.60Mb, 4.06Mb, and 88.84Kb, respectively.

[0039] The simulated whole-genome sequencing data used for ROH detection with different contig N50 reference genome datasets had a depth of 20× and read lengths of 100bp and 150bp, respectively. The BWA mem algorithm is used for comparison, and the ZBOLT SortMarkDup process is used to sort and deduplicate the BAM files. The HaplotypeCaller algorithm in GATK was used to retrieve the original variants for each individual, generating gVCF files. These gVCF files were then merged into a single gVCF file using the CombineGVCFs tool in GATK. Joint variant detection was performed using GenotypeGVCFs in GATK to generate a population-level VCF file. SNPs were extracted from all variants using the parameter "SelectVariants --select-type-to-include SNP" in GATK, followed by hard filtering using "QD<2.0 || MQ<40.0 || FS>60.0 || MQRankSum<-12.5 || ReadPosRankSum<-8.0". The single nucleotide polymorphism (SNP) dataset was further filtered using VCFtools to remove all non-biallelic SNPs and SNPs with a deletion rate exceeding 20%. ROH detection was performed using PLINK software. The VCF file was converted to PLINK file format, and the linkage imbalance (LD) was filtered using the "--indep-pairwise" program, with the window size limited to 50Kb. It is 0.9; The parameters for ROH testing using PLINK software are "--homozyg --homozyg-window-snp 20 --homozyg-density 50 --homozyg-kb 10"; Finally, use F. ROH Error rate assessment at different sequencing depths F ROH The accuracy of the test specifically includes: Calculate the average and total length of the ROH obtained using different sets of simulated data.

[0040] Calculate F using the total length of the ROH. ROH Value, further calculate F ROH Error rate.

[0041] F ROH Error rate limited to .

[0042] ROH detection was performed using reference genomes with different genetic distances, including... The use of real genome data from different species within the same genus, different subspecies within the same species, and different genera within the same family is restricted to be used as reference genomes. Real whole genome data from different subspecies within the same species are used relative to different reference genomes for ROH detection.

[0043] The ROH detection procedure is the same as above.

[0044] This implementation improves the accuracy of ROH detection in the case of low-depth WGS data or low-quality reference genomes by optimizing parameters in PLINK, including: The command "plink --file filename --indep-pairwise n1 n2 n3" is used to calculate the LD between SNP pairs within an n1 KB window, and the process is performed by advancing n2 SNPs in each iteration; if any SNP pair within the window... If the value exceeds n3, one of the SNPs is marked as redundant and deleted; this invention continuously adjusts the window size (n1) of "--indep-pairwise" from 1 Mb to 5 Kb. The threshold (n3) was continuously adjusted from 0.01 to 0.99; The filtered SNP dataset was input into the PLINK software, and ROH detection was performed using the parameters "--homozyg --homozyg-window-snp n4 --homozyg-density n5 --homozyg-kb n6 --homozyg-snp n7". This step uses a sliding window of n4 SNPs to scan the genome, allowing a maximum of n5 SNPs per Mb region to control SNP density, and regions longer than n6 Kb and containing at least n7 SNPs are identified as ROHs. In this invention, the "--homozyg-window-snp" parameter (n4) was continuously adjusted from 100 SNPs to 5 SNPs, and the "--homozyg-snp" parameter (n7) was continuously adjusted from 50 SNPs to 300 SNPs.

[0045] The method proposed in this embodiment is the first to systematically evaluate the impact of key factors such as the contig N50 of the reference genome, the sequencing depth of WGS data, and fragment length on ROH detection, while quantifying the limitations of the genetic distance between the reference genome and the target species on the accuracy of ROH detection.

[0046] This method provides simple and easy-to-use parameter schemes to meet different needs, such as high-quality samples, low-quality samples (e.g., fecal samples), and highly inbred individuals. In particular, it allows for parameter adjustments for situations where data quality may be poor due to the difficulty in sampling endangered species. This helps researchers design sequencing plans in advance based on actual conditions, avoid repeated sequencing, and reduce experimental and time costs.

[0047] Implementation Method 2: The guidance system for detecting animal inbreeding levels based on whole-genome sequencing described in this implementation method includes: The preprocessing unit is used to obtain real high-quality reference genome data and high-quality whole genome sequencing data of a real population from a public database, and to preprocess the obtained data. The preprocessing includes screening complete autosomal sequences with no base deletions from telomeres to telomeres from the high-quality reference genome to construct a new reference genome, and using the new genome as a reference, performing variant detection based on high-quality resequencing data to identify the genetic variations in the genome of each individual in the real population. The haplotype genome file generation unit is used to calculate the inbreeding evaluation parameters ROH and F for different individuals in a real population. ROH The variation data of each individual will be restored to the reference genome, and two haplotype genome files with real variation characteristics will be generated for each individual. The simulation unit is used to simulate a reference genome dataset and a whole genome sequencing dataset using preprocessed real data. The simulated whole genome sequencing dataset includes a dataset with sequencing depth gradient distribution, a dataset with sequencing read length gradient distribution, and a dataset with reference genome continuous gradient distribution. The homozygous fragment detection unit is used to perform homozygous fragment detection using simulated reference genome datasets and whole genome sequencing datasets, respectively, to assess the impact of data quality parameters on detection accuracy. The optimization unit is used to optimize the parameters of the ROH fragment detection software based on the evaluation results.

[0048] Implementation Method 3: A computer device according to this implementation method includes a memory and a processor. The memory stores a computer program. When the processor runs the computer program stored in the memory, the processor executes the guidance method for detecting animal inbreeding levels based on whole-genome sequencing as described in Implementation Method 1.

[0049] Implementation Method 4: A computer-readable storage medium according to this implementation method stores a computer program that, when executed by a processor, performs the steps of the guidance method for detecting animal inbreeding levels based on whole-genome sequencing as described in any one of Implementation Method 1.

[0050] Implementation Method 5, see below Figures 2 to 5 This embodiment describes a specific example of the guidance method for detecting animal inbreeding levels based on whole-genome sequencing as described in Embodiment 1, including: This embodiment aims to address the issue of unclear requirements for sequencing depth and fragment length in ROH detection using WGS data. It simulates WGS data with different sequencing depths and fragment lengths by downloading real Siberian tiger WGS data and constructing haplotype genomes for each WGS dataset for subsequent analysis, and calculates the true F...ROH Value. The specific implementation method is as follows: Step 1: Download WGS data of 13 Siberian tigers and high-quality Siberian tiger reference genome data (GCF_018350195.1, telomere-to-telomere (T2T) level, contig N50=93.60 Mb) from public databases.

[0051] Step 2: From the Siberian tiger reference genome downloaded in Step 1, four T2T autosomes with the highest integrity were selected: A1 (NC_056660.1, approximately 155.91 Mb), B4 (NC_056666.1, approximately 130.82 Mb), D4 (NC_056672.1, approximately 110.53 Mb), and E2 (NC_056674.1, approximately 98.74 Mb). The selected chromosomes contained no deleted base sequences. Using SeqKit, the four autosomes were integrated into a single FASTA format file named "PT-4CHR.fasta," which served as the input file for subsequent simulations of reference genomes with different continuities. The quality of PT-4CHR was evaluated using QUAST, showing a contig N50 of 93.60 Mb, a longest contig length of 155.91 Mb, and no deleted base sequences, meeting the standards for a high-quality reference genome.

[0052] Step 3: Using BWA's mem algorithm, compare the FASTQ format WGS data of 13 Siberian tigers with PT-4CHR to generate a comparison result file in SAM format.

[0053] Step 4: Variation Detection: Using the SAMtools view command, convert the SAM file to BAM format. Call the ZBOLT SortMarkDup workflow to sort the BAM files by chromosome position. Use the GATK HaplotypeCaller algorithm to perform single-sample variation detection on each deduplicated BAM file, generating a gVCF file containing all potential variant sites. Use the GATK CombineGVCFs tool to merge the gVCF files of the 13 samples into a single population-level gVCF file. Then use the GenotypeGVCFs tool to perform population-wide joint variation detection, generating a VCF file containing genotype information for all samples.

[0054] Using GATK's SelectVariants tool, only SNP sites were extracted from the VCF file with the parameter "--select-type-to-include SNP". Then, hard filtering was performed with the filtering criteria set to "QD<2.0 || MQ<40.0 || FS>60.0 || MQRankSum<-12.5 || ReadPosRankSum<-8.0" to remove low-quality SNPs.

[0055] The SNP dataset after hard filtering was further filtered using VCFtools to remove non-bicelesteal SNPs and SNPs with a deletion rate of more than 20%, ultimately obtaining a high-quality population SNP dataset for subsequent ROH detection.

[0056] Step 5: ROH Detection: Convert the high-quality SNP dataset (VCF format) from Step 4 to PLINK binary format using PLINK; perform LD filtering using the "--indep-pairwise" option, with parameters set to a window size of 50Kb and a step size of 25Kb SNPs. The threshold is 0.9. Redundant SNPs with high LD within the window are deleted to avoid interference from LD on ROH detection.

[0057] ROH detection was performed using PLINK's "--homozyg" series of parameters, employing a widely used parameter combination: "--homozyg --homozyg-window-snp 20 --homozyg-density 50 --homozyg-kb 10", which means 20 SNPs within the window, a maximum SNP density of 50 per Mb, and a minimum ROH fragment length of 100 kb, ensuring that the detected ROH fragments are biologically significant. (The formula is used.) The inbreeding coefficient was calculated for each sample, and only ROH fragments longer than 100 kb were counted, excluding short fragments of random homozygous pairs F. ROH The impact.

[0058] Step 6: Select 8 tigers (F) from 13 Siberian tigers. ROH Individuals whose values ​​exhibit a clear gradient distribution (F) ROH = 7.04%, 13.98%, 17.91%, 23.28%, 29.13%, 36.43%, 38.84%, 48.36%), which will be used as the actual template for subsequent WGS data simulation. The selection results are shown in the table below: Table 1F ROH Real data of individuals with value gradient distribution

[0059] Step 7: Using the PT-4CHR genome as the input file, and based on the VCF files of the 8 individuals generated in Step 4 above, replace alleles using a custom script to generate two haplotype genome files for each individual. The specific rules are as follows: For the "1 / 1" genotype (homozygous variant, different from PT-4CHR): alleles in PT-4CHR are replaced with alternative alleles in VCF, and both haplotype genomes are replaced; for the "0 / 1" genotype (heterozygous variant): allele replacement is performed only in one of the two haplotype genomes; for the "0 / 0" genotype (homozygous congruent): no replacement is performed, maintaining congruence with PT-4CHR. Finally, 16 FASTA format haplotype genomes are obtained, named "ATxx-hap1", "ATxx-hap2", etc.

[0060] Step 8: WGS data simulation with different sequencing depths and fragment lengths Using wgsim software, we simulated basic WGS data with 16 haplotype genomes as input files, with parameters set to "-1 100 -2 100 -e 0" (-1 / -2 represents paired-end read length, and -e 0 ignores sequencing errors). Each haplotype genome generated approximately 100 Gb of data (equivalent to 200× coverage of the PT-4CHR genome) for subsequent analysis.

[0061] The baseline WGS data was randomly sampled using SeqKit to generate datasets at 15 sequencing depths: 2×, 2.5×, 3×, 3.5×, 4×, 5×, 6×, 8×, 10×, 12×, 15×, 20×, 30×, 60×, and 100×. Each depth corresponds to paired-end FASTQ files for 8 individuals, with a read length limit of 100 bp. Simulated data at depth 1× was excluded because the ROH fragment could not be detected.

[0062] WGS data with different read lengths were simulated using wgsim. The read lengths were set to 30bp, 50bp, 75bp, 100bp, and 150bp (shorter read lengths represent ancient DNA or samples of lower quality; 100bp and 150bp are commonly used read lengths). The parameters were adjusted to "-1 read length -2 read length -e 0". Each read length corresponds to paired-end FASTQ files of 8 individuals. The final depth was limited to 20×.

[0063] Step 9: Using PT-4CHR as a reference genome, ROH detection was performed using WGS simulated data at different depths with read lengths of 100 bp. ROH was detected using standard PLINK parameters. Results showed that only when the individual sequencing depth reached or exceeded 15-fold did the detectable ROH fragments in the genome closely approximate the true state. Figure 3 As shown. (Through F) ROH Error rate to test F ROH Detection accuracy:

[0064] At a sequencing depth of 15x, F ROH The error rate was only 0.84% ​​± 1.71%, then gradually decreased to 0.03% ± 0.08% at 100x, and lower depths may underestimate (<3×) or overestimate (3×~15×) F. ROH This deviation in F ROH This was more pronounced in <30% of individuals. However, when F ROH When the inbreeding rate is >30%, a sequencing depth of 3 times or more can accurately reflect the true level of inbreeding, such as... Figure 4 As shown. Therefore, the recommended minimum sequencing depth is 15×, and even for samples with poor DNA quality, a sequencing depth of 5× should be achieved, but F ROH It may be overestimated.

[0065] Step 10: Based on the results of Step 8, perform ROH detection on data of different fragment lengths using 20×WGS simulated data, and calculate F. ROH When using the same reference genome, the F-simulation read lengths of 100bp and 150bp are different. ROH There was no significant difference between the average ROH length and the data, indicating that a sequencing read length of 100 bp is generally sufficient for ROH detection. At shorter read lengths (30 bp, 50 bp, and 75 bp), although ROH detection improved with increasing read length, the effect of read length on ROH detection was not significant (e.g., ...). Figure 5 (As shown), however, ROH detection results in ancient DNA studies still need to be interpreted with caution because: sequencing depth is usually very low in ancient DNA studies; DNA damage exists in ancient DNA sequences; and many reads shorter than 30 bp exist in ancient materials.

[0066] Implementation Method Six: This implementation method provides a specific example of the guidance method for detecting animal inbreeding levels based on whole-genome sequencing as described in Implementation Method One, including: This embodiment aims to address the issue of unclear requirements for reference genome continuity in ROH detection. It constructs a new high-quality reference genome by downloading a high-quality reference genome and screening chromosomes without base deletions. Deletion sequences are then inserted to simulate different continuity reference genomes, and ROH detection is performed using reference genomes with different contig N50 values. The specific implementation method is as follows: Step 1: Using a custom script, 300 bp consecutive "N" sequences are randomly inserted as deletion sequences into the PT-4CHR genome constructed in Implementation Method 5 to simulate base deletions in real genome assembly. The core logic of the script is: on each chromosome of the PT-4CHR, the deletion sequence is inserted according to the principle of "random position and non-overlapping". After each insertion, the Contig N50 is calculated until a reference genome with the target gradient continuity is obtained. Finally, 13 simulated reference genomes with different Contig N50s are generated, ranging from 25.91 Kb to 93.60 Mb. The specific parameters are shown in the table below: Table 2 shows simulated reference genomes with different contig N50 values.

[0067] Step 2: ROH detection was performed using simulated reference genomes with different continuities. The detection method and parameter settings were consistent with Example 1, using simulated WGS data of 20×100bp as input and PLINK software for ROH detection. When the contig N50 of the reference genome decreased from 93.60Mb to 1.2Mb, the total number of detectable ROH fragments remained stable; when the contig N50 decreased to below 1.2Mb, the number of detectable ROH fragments increased significantly. This indicates that as genome continuity decreases, larger ROH fragments are more easily fragmented during detection. Figure 2 As shown. When contig N50 reaches 4Mb, F can be accurately evaluated. ROH The average length of ROH fragments decreased as contig N50 decreased, as shown in Table 3.

[0068] Table 3. Average ROH length for ROH detection using different continuous reference genomes.

[0069] (Continued from Table 3)

[0070] Implementation Method Seven, see below Figure 6 and Figure 7This embodiment describes a specific example of the guidance method for detecting animal inbreeding levels based on whole-genome sequencing as described in Embodiment 1, including: This embodiment aims to address the unclear impact of genetic distance between the reference genome and the target population on ROH detection. By selecting genomes of closely related species at different genetic distances as references, ROH detection is performed on real WGS data of three tiger subspecies to assess the influence of genetic distance on ROH detection. The specific implementation method is as follows: Step 1: Download high-quality reference genomes of five feline species from the NCBI / CNGB database, following the genetic distance gradients of subspecies, genus, and family. Specific information is shown in the table below: Table 4 Reference Genome Information for Different Felids

[0071] Step 2: Download data from a public database of 40 tiger groups, including the 13 Siberian tigers used in Implementation Method 5 (also including 12 Bengal tigers and 15 South China tigers). The depth is between 14.55× and 27.57×, the coverage is between 98.39% and 99.87%, and the data quality meets the analysis requirements.

[0072] Step 3: Using the method described in Implementation Method 5, WGS data from three tiger subspecies populations were mapped onto high-quality reference genomes of five feline species. After variation detection and data filtering, ROH detection was performed.

[0073] Step 4: Compare the differences in ROH test results using different indicators: F ROH Value according to formula: Calculate, where ∑ L ROH It is the total length of the detected ROH fragments. Lauto This represents the total length of autosomes in the corresponding reference genome; The percentage of ROH fragments of different lengths, namely, the percentages of ROH fragments >25Mb, 10~25Mb, 5~10Mb, 1~5Mb, 500kb~1Mb, and less than 500kb.

[0074] The results showed that the F-values ​​of the three tiger subspecies populations were calculated using reference genomes from the three tiger subspecies. ROH Highly consistent, the F-strain of the three tiger subspecies populations was calculated using the lion and domestic cat genomes as reference genomes. ROHThe number of ROH fragments detected was significantly lower than that calculated based on three tiger subspecies. Length distribution analysis of the detected ROH fragments showed significant differences in ROH detection results below the 5Mb threshold. For ROH fragments smaller than 5Mb, the number of detectable ROH fragments increased with increasing genetic distance. The number of ROH fragments larger than 5Mb detected using the three tiger subspecies reference genomes was significantly higher than the number detected using the lion and domestic cat genomes. For ROH fragments longer than 5Mb, using the Bengal tiger genome as a reference, the number of ROH fragments detected in the Siberian tiger and South China tiger populations was lower than that in the Siberian tiger and South China tiger genomes. Figure 6 As shown, this indicates that reference genomes within subspecies should be prioritized in the analysis, but the influence of genetic distance cannot be completely avoided.

[0075] Based on the above three embodiments, 5 to 7, the present invention forms a practical method for ROH detection using whole-genome sequencing technology in protective genomics, such as... Figure 7 As shown.

[0076] Implementation Method Eight: This implementation method provides a specific example of the guidance method for detecting animal inbreeding levels based on whole-genome sequencing as described in Implementation Method One, including: This embodiment aims to address the problem that the default PLINK parameters are not suitable for low-quality data. By setting different combinations of PLINK parameters, ROH detection is performed on WGS data with low depth (5×) and low continuity reference genome (contig N50=0.09 Mb), and the optimal parameter combination is screened to improve the detection accuracy of low-quality data.

[0077] Step 1: Select WGS simulated data from 8 Siberian tigers in Implementation Method 5 with a sequencing depth of 5× and a read length of 100 bp to represent low-quality samples such as feces and ancient DNA; select the PT-4CHR-11 genome (contig N50=0.09 Mb) simulated in Implementation Method 6 to represent the fragmented reference genome.

[0078] Step 2: Focusing on the two core steps of PLINK detection of ROH (SNP filtering and ROH scanning), adjust the parameters for the low-depth WGS data in Implementation Method 5. The parameter adjustment also follows the principle of gradient setting.

[0079] When the parameter "--homozyg-window-snp" is reduced from 100 to 20, F ROH Both the average ROH length and the length of the ROH will increase; looking at them separately, for F ROH Calculations show that setting this parameter to 80 will allow most individuals to achieve F... ROH The value is close to the true value, while 70 is more suitable for highly inbred individuals (F). ROH(≈50%). For calculating the average ROH length, it is recommended to use "--homozyg-window-snp20", but this may underestimate the ROH length.

[0080] The parameters "--homozyg-snp 150" and "--indep-pairwise 50 1 0.8" can also improve the accuracy of ROH detection at low sequencing depths.

[0081] For low-sequencing-depth WGS data, this embodiment recommends using "--indep-pairwise 50 1 0.8; --homozyg-window-snp 80; --homozyg-snp 150". For F... ROH For more than 50% of individuals, it would be more appropriate to reduce "--homozyg-window-snp" to 70.

[0082] Step 3: Parameter adjustment for the low continuity reference genome in Implementation Method 6, mainly targeting low-quality reference genomes with contig N50 < 100kb.

[0083] Adjusting the "--homozyg-window-snp" parameter reduces the number of SNPs in the window, as some ROH fragments may span multiple contigs, leading to false negatives. Lowering this parameter increases the number of detectable ROH fragments, but has no significant effect on the average ROH length. For reference genomes with contig N50 < 100 Kb, this example recommends setting "--homozyg-window-snp" to 10-15, while keeping other parameters at optimal values, such as "--homozyg-snp 150" and "--indep-pairwise 50 1 0.8".

[0084] Step 4: To address the adverse effects of the large genetic distance between the reference genome and the target population described in Implementation Method 7, parameters were adjusted. The ROH values ​​in the Bengal tiger population WGS data were analyzed using the lion reference genome. It was found that smaller values ​​for "--indep-pairwise" and "--homozyg-window-snp" resulted in lower F... ROH The average ROH length increases, but is positively correlated with "--indep-pairwise" and negatively correlated with "--homozyg-window-snp". Due to taxonomic unit-specific differences, this embodiment does not recommend using specific parameters, but for large felines, it is suggested to evaluate F. ROHWhen evaluating ROH length, the values ​​of "--indep-pairwise" and "--homozyg-window-snp" should be decreased, while the value of "--indep-pairwise" should be increased.

[0085] Those skilled in the art will understand that embodiments of this disclosure can be provided as methods, systems, or computer program products. Therefore, this disclosure can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this disclosure 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, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0086] This disclosure is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this disclosure. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create a machine for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 The computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to operate in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes. These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0087] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this disclosure and not to limit its protection scope. Although this disclosure has been described in detail with reference to the above embodiments, those skilled in the art should understand that after reading this disclosure, they can still make various changes, modifications or equivalent substitutions to the specific implementation of the invention, but these changes, modifications or equivalent substitutions are all within the protection scope of the published pending claims.

Claims

1. A guide method for detecting inbred level of an animal based on whole genome sequencing method, characterized in that, The method comprises: Step S1: obtaining real high-quality reference genome data and high-quality whole genome sequencing data of a real population from a public database, and preprocessing the obtained data, the preprocessing comprising constructing a new reference genome from the high-quality reference genome by screening complete autosomal sequences with no base deletion characteristics from telomere to telomere, and performing variation detection based on the high-quality resequencing data to identify genetic variations of each individual in the real population on the genome by taking the new genome as a reference; Step S2: Calculate inbreeding evaluation parameters ROH and F of different individuals in the real population ROH and restore the variant data of each individual to the reference genome, and each individual will generate 2 haplotype genome files with true variant characteristics; Step S3: using the preprocessed real data to simulate reference genome data sets and whole genome sequencing data sets, wherein the simulated whole genome sequencing data sets comprise data sets with gradient distribution of sequencing depth, data sets with gradient distribution of sequencing read length, and data sets with gradient distribution of reference genome continuity; Step S4: performing homozygosity segment detection using the simulated reference genome data sets and whole genome sequencing data sets respectively to evaluate the influence of data quality parameters on detection accuracy; Step S5: optimizing parameters of ROH segment detection software under low-quality sequencing data based on the evaluation results.

2. The method according to claim 1, wherein the method is characterized by, The F ROH The calculation method is: Among them, F ROH This represents the ratio of the total physical length of autosomes. This represents the total length of the ROH segment on an autosome. The total length of the autosome is 100 kb; the ROH fragment length threshold is set to 100 kb.

3. The method according to claim 1, wherein the method is characterized by, The sequencing depth gradient comprises 100x, 60x, 30x, 20x, 15x, 12x, 10x, 8x, 6x, 5x, 4x, 3.5x, 3x, 2.5x, 2x and 1x; the sequencing read length gradient comprises 30bp, 50bp, 75bp, 100bp and 150bp; and the reference genome is characterized by the value of contig N50, and the gradient in the reference genome comprises 25.91Kb, 44.90Kb, 88.84Kb, 177.75Kb, 304.10Kb, 600.66Kb, 1.20Mb, 2.05Mb, 4.06Mb, 8.13Mb, 13.89Mb, 27.38Mb, 54.88Mb and 93.60Mb.

4. The method according to claim 1, wherein the method is characterized by, In step S3, the data in the simulated whole genome sequencing data set is generated using the wgsim software, and the parameter setting is fixed read length without introducing sequencing errors; and the haplotype genome file is generated by replacing the variant sites in the reference genome, wherein when the variant site is 1 / 1 genotype, the bases of the two haplotypes at this position are replaced with the allele, and when the genotype is 0 / 1, only one haplotype is replaced.

5. The method according to claim 1, wherein the method is characterized by, The homozygosity segment detection in step S4 comprises: Using BWA software mem An algorithm aligns whole genome sequencing data to reference gene sequencing data; The generated BAM file is preprocessed, including sorting and deduplication processing; The preprocessed BAM file is subjected to population-level variation detection, including; The HaplotypeCaller algorithm in the GATK toolkit is used to call original variations of each individual, and a gVCF file is generated; The gVCF files are combined into one gVCF file using the CombineGVCFs tool in the GATK toolkit, joint variation detection is performed using GenotypeGVCFs, and a population-level VCF file is generated; Single nucleotide polymorphisms (SNPs) are extracted from all variations, and the single nucleotide polymorphisms (SNPs) are hard filtered; and the single nucleotide polymorphisms (SNPs) are further filtered by VCFtools software to remove all non-biallelic SNPs and SNPs with a deletion rate of more than 20%; ROH detection was performed using PLINK software, and the VCF file was converted into PLINK file format; linkage disequilibrium (LD) filtering was performed using the --indep-pairwise command, and the window size was set to 50Kb, and the correlation index was set to 0.9 a threshold value of 0.9; Homozygous fragment detection is performed using the --homozyg command, including: Each sliding window contains at least 20 SNPs, the SNP density in the window is one out of every 50 SNPs, and the physical length of the ROH fragment is at least 100Kb.

6. The method according to claim 1, wherein the method is characterized by, The step S4 further comprises using F ROH Error rate evaluation of F ROH Accuracy of detection, F ROH The error rate calculation formula is: .

7. The method according to claim 1, wherein the method is characterized by, The step S5 comprises: The command plink --file filename --indep-pairwise n1 n2 n3 is used to calculate the linkage disequilibrium LD between SNP pairs within the window of n1, and the processing is performed in a manner of advancing n2 SNPs per iteration; if the correlation index of any SNP pair within the window is greater than n3, one of the SNPs is marked as redundant and deleted; the window size n1 of the command --indep-pairwise is continuously adjusted from 1 Mb to 5 Kb, and the threshold n3 is continuously adjusted from 0.01 to 0.

99. The command --homozyg --homozyg-window-snp n4 --homozyg-density n5 --homozyg-kb n6 --homozyg-snp n7 is used for ROH detection, which uses a sliding window of n4 SNPs to scan the genome, allows a maximum of n5 SNPs per Mb region to control the SNP density, and identifies a region with a length exceeding n6 and containing at least n7 SNPs as ROH; the parameter n4 of the command --homozyg-window-snp is continuously adjusted from 100 SNPs to 5 SNPs, and the parameter n7 of the command --homozyg-snp is continuously adjusted from 50 SNPs to 300 SNPs. The command --homozyg --homozyg-window-snp n4 --homozyg-density n5 --homozyg-kb n6 --homozyg-snp n7 is used for ROH detection, which uses a sliding window of n4 SNPs to scan the genome, allows a maximum of n5 SNPs per Mb region to control the SNP density, and identifies a region with a length exceeding n6 and containing at least n7 SNPs as ROH; the parameter n4 of the command --homozyg-window-snp is continuously adjusted from 100 SNPs to 5 SNPs, and the parameter n7 of the command --homozyg-snp is continuously adjusted from 50 SNPs to 300 SNPs.

8. A guide system for detecting inbred level of an animal based on whole genome sequencing method, characterized in that, The system comprises: A preprocessing unit is configured to obtain real high-quality reference genome data and high-quality whole genome sequencing data of a real population from a public database, and to preprocess the obtained data, the preprocessing comprising constructing a new reference genome from the high-quality reference genome by screening complete autosomal sequences with no base deletion from telomere to telomere, and performing variation detection based on the high-quality resequencing data with the new genome as a reference to identify genetic variations of each individual in the real population on the genome; A haplotype genome file generating unit for calculating inbreeding evaluation parameters ROH and F of different individuals in a real population ROH And restore the variation data of each individual to the reference genome, and each individual will generate 2 haplotype genome files with true variation characteristics; An analog unit is configured to simulate reference genome data sets and whole genome sequencing data sets using the preprocessed real data, the simulated whole genome sequencing data sets including data sets with gradient distribution of sequencing depth, data sets with gradient distribution of sequencing read length, and data sets with gradient distribution of reference genome continuity; A homozygous fragment detection unit is configured to perform homozygous fragment detection using the simulated reference genome data sets and whole genome sequencing data sets respectively to evaluate the influence of data quality parameters on detection accuracy; An optimization unit is configured to optimize parameters of ROH fragment detection software under low-quality sequencing data based on the evaluation results.

9. A computer device, comprising: The computer readable storage medium stores a computer program, and when the processor runs the computer program stored in the memory, the processor executes the guidance method for detecting the inbreeding level of an animal based on whole genome sequencing according to any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer readable storage medium stores a computer program, and when the processor runs the computer program stored in the memory, the processor executes the guidance method for detecting the inbreeding level of an animal based on whole genome sequencing according to any one of claims 1-7.

Images