A method for ancestral haplotype typing across tetraploid and hexaploid wheat at genomic fine intervals

Abstract

The application discloses a method for ancestral haplotype typing in a fine genomic interval across tetraploid and hexaploid wheat. The application provides a method, steps of which are as follows: 1) obtaining whole genome resequencing data of each wheat sample to be tested; 2) obtaining filtered DNA sequence alignment files of each sample to be tested; 3) SNP identification and CNV variation identification; 4) obtaining SNP polymorphism site density distribution; 5) fitting the density distribution of the SNP polymorphism sites to obtain a normal distribution, and determining a genetic distance hard threshold according to a boundary; 6) obtaining an initial ancestral haplotype; 7) obtaining an initial priority order of the sample; and 8) obtaining an ancestral haplotype file. The application does not need resequencing data of a reference population as input, which reduces the total detection cost. The application also does not need to limit the phylogenetic relationship between populations to which the samples to be detected belong, which can improve the detection accuracy in crops with complex evolutionary history.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of plant genomics, and particularly relates to a method for ancestral haplotype typing in a fine genomic interval across tetraploid and hexaploid wheat. BACKGROUND

[0002] The progress of crop breeding largely depends on the utilization of variation in the major gene pool of the species, while a large amount of genetic diversity of major crops is lost during domestication. Therefore, exploring new and underutilized genetic resources is an important part of crop breeding. There are secondary gene pools represented by wild relatives for major crops such as wheat, maize, and rice. In recent years, extensive studies have shown that some introgression fragments from wild species, ancestral species, and different lineages have the potential to enhance crop adaptability.

[0003] Therefore, systematically identifying the ancestral haplotype type of crop genomic fragments and introducing haplotype genomic fragments that have not been extensively utilized in wild species and relatives can greatly expand the available genetic variation resources and is of great significance to crop breeding. Moreover, clarifying the ancestral origin of crop genomic fragments can also promote the study of crop evolution rules. However, the current mainstream method for detecting distant genomic introgression fragments requires a reference population as a control, and has strict requirements for the phylogenetic relationship between the populations to be detected, which limits its application range in crops with complex evolutionary history.

[0004] Common wheat is an important food crop and is an allohexaploid with a large genome and high repetitive sequences. With the publication of the common wheat reference genome and the reduction of sequencing costs based on the second-generation sequencing, large-scale resequencing of varieties, relatives, and wild species has become an important technology and method for systematically analyzing the evolution of germplasm resources. Large-scale resequencing can be used to obtain comprehensive and high-quality genotype data, thereby providing the possibility for accurately inferring the ancestral haplotype type of blocks at the whole-genome level. SUMMARY

[0005] An object of the present application is to provide a method for determining the ancestral haplotype to which a genomic block of different wheat samples belongs across tetraploid and hexaploid.

[0006] The method provided by the present application comprises the following steps:

[0007] 1) obtaining whole-genome resequencing data of each wheat sample to be tested;

[0008] In step 1), the data can be obtained from a database or by whole-genome sequencing.

[0009] The wheat sample is a variety or strain, a relative or wild species;

[0010] The ploidy of the wheat sample is tetraploid or hexaploid.

[0011] 2) Filter the whole genome resequencing data of each wheat sample to be tested to achieve quality control, and obtain filtered data; then align the filtered data with the reference wheat genome sequence to obtain a DNA sequence alignment file of each sample to be tested; filter the DNA sequence alignment file of each sample to be tested again to obtain a filtered DNA sequence alignment file of each sample to be tested;

[0012] Specifically, the software used for the above filtering to achieve quality control is (http: / / www.usadellab.org / cms / index.php?page=trimmomatic), and the specific parameters are LEADING: 3 (reads with a base quality value lower than the threshold value at the beginning are discarded); TRAILING: 3 (reads with a base quality value lower than the threshold value at the end are discarded); SLIDINGWINDOW: 4:15 (sliding window size starting from the 5' end of reads); MINLEN: 36 (length of reads after cutting below the threshold value are discarded), and finally the filtered data is obtained for subsequent analysis;

[0013] The above alignment method is to use the BWA-MEM tool in the BWA software (https: / / sourceforge.net / projects / bio-bwa / files / ) to align the filtered data on the Chinese Spring wheat reference genome (IWGSC RefSeqv1.0, https: / / urgi.versailles.inra.fr / download / iwgsc / IWGSC_RefSeq_Annotations / v1.1 / ), and select reads that exist in the "best hit".

[0014] The above re-filtering is performed using bamtools (v2.5.1, https: / / github.com / pezmaster31 / bamtools) and samtools (v1.4, http: / / www.htslib.org / ) tools: using bamtools to remove reads with a distance between read pairs greater than 10,000 bp or 0 bp and reads with a back-posting quality less than 1, and using samtools to remove PCR duplicate data.

[0015] Through the above alignment and re-filtering, the BAM file is finally obtained for subsequent analysis, denoted as the filtered DNA sequence alignment file of each sample to be tested.

[0016] 3) Perform SNP calling and CNV (copy number variation) calling on the filtered DNA sequence alignment files of all the samples to be tested, to obtain a filtered VCF file recording all SNP site information and a BED file containing CNV variation interval information for each sample to be tested;

[0017] The software for SNP calling is GATK software, and the SNP site filtering parameters are set to "QD < 2.0, FS > 60.0, MQRankSum < -12.5, ReadPosRankSum < -8.0, SOR > 3.0, MQ < 40.0, DP > 30 || DP < 3".

[0018] The CNV variation calling is to divide the reference whole genome into small windows with 5Mbp as a unit, and use bedtools (v2.26.0, https: / / bedtools.readthedocs.io / en / latest / ) software to calculate the "average coverage depth" (Dep bin ) of the resequencing alignment reads (i.e. the DNA sequence alignment file of each sample to be tested) in each window of each sample to be tested; ave ) of the whole genome of the material (calculated from the DNA sequence alignment file of each sample to be tested); and normalize by dividing the two to obtain the "average relative coverage depth" of each small window: Dep bin / Dep ave . Dep bin / Dep ave Windows with a value lower than 0.5 or higher than 1.5 are considered "CNV variation intervals", and the BED file containing CNV variation interval information for each sample to be tested is finally obtained.

[0019] 4) Find and count the number of different SNPs in the small windows that are not CNV variation intervals for each of the two samples in the filtered VCF file recording all SNP site information and the BED file containing CNV variation interval information for each sample to be tested, and divide the number of each SNP by the length of the corresponding small window to obtain the SNP density, to obtain a genetic distance matrix file recording the SNP density between the two varieties in the small window; and plot the SNP density of all samples to obtain the distribution of SNP polymorphic site density;

[0020] The method 4) of plotting the SNP density of all samples to obtain the distribution of SNP polymorphic site density comprises the following steps:

[0021] 4-1) Select a fixed interval length as a unit to divide the reference genome into small windows within 1-5Mbp;

[0022] 4-2) Using PLINK v1.9 software, count the number of different SNPs in each pair of samples in the small windows where neither sample is in the CNV variant interval, divide each SNP number by the length of the corresponding small window to get the SNP density, and obtain the genetic distance matrix file recording the SNP density between the two varieties in the small window;

[0023] 4-3) Plot the SNP density of all samples by taking the logarithm with base 10, and obtain the density distribution.

[0024] 5) Fit the density distribution of the SNP polymorphic sites to obtain a normal distribution, and determine the genetic distance hard threshold according to the boundary;

[0025] In step 5), the fitting algorithm is the maximum likelihood algorithm.

[0026] The hard threshold is determined by taking the cross boundary of the fitted normal distribution as the hard threshold.

[0027] In the embodiments of the present application, for wheat, the hard threshold is 10 -3 per bp.

[0028] 6) Cluster the genetic distance matrix of each small window obtained in step 4) to obtain a hierarchical clustering tree sample for each small window, and use the genetic distance hard threshold obtained in step 5) to divide the hierarchical clustering tree sample of each small window by height to obtain an initial ancestral haplotype;

[0029] 7) Calculate the initial ancestral haplotype to obtain a sample subset that maximizes the diversity of ancestral haplotypes, and then add the samples not selected into the subset in random order to the subset; the order in which the samples are selected into the subset is used as the initial priority order of the samples;

[0030] 8) Calculate the initial ancestral haplotype obtained in step 6) and the sample priority order obtained in step 7) to obtain a sample priority order file and an ancestral haplotype file; the purpose of the calculation is to minimize the total number of changes in ancestral haplotypes between adjacent small windows.

[0031] In the above step 8), the calculation method can be any one of the following two modes:

[0032] Unsupervised mode, i.e., the priority order of the samples is completely determined by the class bubble algorithm;

[0033] Semi-supervised mode, i.e., the priority order of the samples is determined by the class bubble algorithm within a pre-defined sample grouping.

[0034] In the embodiments of the present application, a semi-supervised mode is adopted.

[0035] Specifically:

[0036] The calculation (sample priority ordering method) in the above method step 8) is either of the following two modes:

[0037] 8-1) Unsupervised mode, the priority order of the samples is changed using a bubble sort algorithm so that the total number of ancestral haplotype changes between adjacent small windows is minimized, starting from the sample with the highest priority, two samples with adjacent priority are exchanged, and the ancestral haplotype is determined according to the priority order of the samples before and after the exchange, respectively, and the total number of ancestral haplotype changes is calculated, if the total number of changes is reduced, the exchange is retained, otherwise the priority order is restored. Until the exchange of sample order cannot reduce the total number of haplotype changes, the step is ended. Then the ancestral haplotype type of all small windows is determined according to this sample priority order to obtain the ancestral haplotype file.

[0038] After determining the sample priority order in each round, the ancestral haplotype determination method is as follows 1)-3):

[0039] 8-1-1) Assign the same new ancestral haplotype type to all small windows of the sample with the highest priority;

[0040] 8-1-2) Starting from the sample with the second highest priority, for each small window of each sample, if there is a sample with a higher priority in the same initial ancestral haplotype as the sample in the small window, the ancestral haplotype type of the sample in the small window is determined as the same ancestral haplotype type as the high-priority sample;

[0041] 8-1-3) Starting from the sample with the second highest priority, for each small window of each sample, if there is no sample with a higher priority in the same initial ancestral haplotype as the sample in the small window, the ancestral haplotype type of the sample in the small window is determined as a new ancestral haplotype type.

[0042] After determining the ancestral haplotype type, the calculation method of the total number of ancestral haplotype changes is: for each sample, if the ancestral haplotype types of adjacent small windows are different, the total number of ancestral haplotype changes is increased by 1.

[0043] 8-2) Semi-supervised mode, each sample should belong to a pre-defined group. Within each group, the priority order of the samples is determined by the unsupervised mode. Then the priority order of the samples is combined according to the priority order between the groups to obtain the total priority order, and the ancestral haplotype type of all small windows is determined according to the sample priority order to obtain the ancestral haplotype file.

[0044] The method further comprises the following step 9): correcting the ancestral haplotype file to obtain a corrected ancestral haplotype file.

[0045] In the method, in step 9), the correction is performed by a Bayesian method to fine-tune the interval division and reduce random errors caused by a hard threshold.

[0046] The method further comprises the following step 10): performing visual processing on the corrected ancestral haplotype file to obtain a visualized map of the ancestral haplotype attribution state of each segment on the whole genome of each sample to be tested; and determining, according to segments with the same color in different samples to be tested in the visualized map, that the segment is an ancestral haplotype block.

[0047] The ancestral haplotype is a group of haplotypes with similar single nucleotide polymorphisms inherited from a common ancestor.

[0048] The ancestral haplotype block is a genomic interval with the same ancestral haplotype type that is not broken by crossover recombination.

[0049] The technical problem to be solved by the present application is how to infer the ancestral haplotype to which a wheat genomic block belongs. The main flowchart of the present application is as follows Figure 1 The present application does not need resequencing data of a reference population as input, which reduces the total cost of detection. The present application also does not need to limit the phylogenetic relationship between the population to which the sample to be detected belongs, which can improve the detection accuracy in crops with complex evolutionary history. BRIEF DESCRIPTION OF DRAWINGS

[0050] Figure 1 The figure is a flowchart of the method.

[0051] Figure 2 The figure is a SNP density distribution map of two samples in a small window of the whole genome.

[0052] Figure 3 The figure is two normal distributions calculated using an EM algorithm and a threshold for distinguishing between the two.

[0053] Figure 4 The figure is a whole genome visualization result of the ancestral haplotype distribution of four wheat materials.

[0054] Figure 5 The figure is a PCR product sequencing verification result of the interval with the same ancestral haplotype type of Lao Laihe and Hong Chunmai and the different ancestral haplotype type of Loflin 10. DETAILED DESCRIPTION

[0055] In the following examples, the experimental methods used are conventional methods unless otherwise specified.

[0056] The materials, reagents, etc. used in the following examples, unless otherwise specified, can be obtained commercially.

[0057] Example 1, Establishment of a method for inferring the ancestral haplotypes to which the blocks of the wheat genome belong

[0058] 1. Obtaining whole genome resequencing data of wheat

[0059] The whole genome resequencing data of hexaploid wheat Laoxiala, Hongchunmai, and Lovrin 10 were downloaded from the National Center for Biotechnology Information database (NCBI) using index number PRJNA596843, and the whole genome resequencing data of tetraploid wheat Zavitan were downloaded using index number PRJEB31422, i.e. the original sequencing data of each sample to be tested (i.e. the whole genome resequencing data of each wheat sample to be tested).

[0060] 2. Quality control of resequencing data and genome alignment

[0061] 1) Quality control of resequencing data

[0062] The original sequencing data of each sample to be tested obtained in 1 above was filtered to achieve quality control using Trimmomatic v0.36 software, the software source for filtering was (http: / / www.usadellab.org / cms / index.php?page=trimmomatic), low-quality sequencing data was removed, and the specific parameters were LEADING: 3 (reads starting base quality value below the threshold value was discarded); TRAILING: 3 (reads ending base quality value below the threshold value was discarded); SLIDINGWINDOW: 4:15 (sliding window size from the 5' end of reads); MINLEN: 36 (length of reads after cutting below the threshold value was discarded), and finally the filtered data was used for subsequent analysis.

[0063] 2) Genome alignment of resequencing data

[0064] The filtered data was aligned on the Chinese Spring wheat reference genome (IWGSCv1 RefSeq v1.0, https: / / urgi.versailles.inra.fr / download / iwgsc / IWGSC_RefSeq_Annotations / v1.1 / ) using the BWA-MEM tool in the BWA software (https: / / sourceforge.net / projects / bio-bwa / files / ), and read pairs with the "best hit" were selected for retention, to obtain the DNA sequence alignment file (BAM file) of each sample to be tested.

[0065] The BAM files were then re-filtered using bamtools (v2.5.1, https: / / github.com / pezmaster31 / bamtools) and samtools (v1.4, http: / / www.htslib.org / ) tools. Bamtools was used to remove reads with a distance between read pairs greater than 10,000 bp or 0 bp and reads with a mapping quality less than 1. Samtools was used to remove PCR duplicate data. The re-filtered BAM files were finally obtained for subsequent analysis, denoted as the filtered DNA sequence alignment files of each test sample.

[0066] 3. Filtering and identifying variant sites and intervals of the test wheat materials

[0067] 1) Filtering and identifying SNP variant sites

[0068] The functions of HaplotypeCaller, GenotypeGVCF, SelectVariants and VariantFiltration of GATK software (v3.8, https: / / gatk.broadinstitute.org / hc / en-us) were used to identify single nucleotide polymorphism (SNP) of all the filtered DNA sequence alignment files of the test samples in the above 2, and the filtered VCF files recording all the SNP site information were obtained.

[0069] In the above single nucleotide polymorphism (SNP) identification, the SNP site filtering parameters are set as "QD < 2.0, FS > 60.0, MQRankSum < -12.5, ReadPosRankSum < -8.0, SOR > 3.0, MQ < 40.0, DP > 30 || DP < 3". The meanings of each parameter are as follows: QD: QualByDepth is the quality score of the allele site depth normalized; FS: FisherStrand chain direction Fisher's exact test value; MQRankSum: MappingQualityRankSumTest back-post quality rank sum test value; SOR: StrandOddsRatio chain odds ratio; ReadPosRankSum: Read Position Rank Sum Test read position rank sum test; MQ: RMSMappingQuality of all read segments of the site back-post quality root mean square; DP: Depth of the site read depth.

[0070] 2) Filter and identify CNV variant segments

[0071] The reference genome is divided into small windows with 5Mbp as a unit, and the "average coverage depth" (Dep bin ) of the filtered DNA sequence alignment file of each sample obtained in the above 2 is calculated in each window using bedtools (v2.26.0, https: / / bedtools.readthedocs.io / en / latest / ) software; and the "average read coverage depth" (Dep ave ) of the whole genome of the material is calculated by the "average coverage depth" of the material in all small windows of the whole genome; and the two are divided to obtain the "average relative coverage depth" of each small window: Dep bin / Dep ave . Dep bin / Dep ave The small window with a value lower than 0.5 or higher than 1.5 is regarded as a "CNV variant interval", and finally the BED file containing CNV variant interval information of each sample is obtained.

[0072] 4, Calculate SNP variant site density distribution

[0073] 1) First, divide the reference genome into small windows with 5Mbp as a unit;

[0074] 2) Using PLINK v1.9 software, from the filtered VCF file containing all SNP locus information obtained in 3.1) above and the BED file containing CNV variation interval information for each sample to be tested obtained in 3.2) above, the number of differential SNPs in the small window where neither sample is a CNV variation interval is counted for each pair of samples (any two varieties are combined into one pair of samples). The density of each SNP is obtained by dividing the number of each SNP by the corresponding small window length, and a genetic distance matrix file containing the density of differential SNPs between the two varieties in the small window is obtained.

[0075] 3) Based on the logarithm of the SNP density of all samples, plot the density distribution using R software (https: / / cran.r-project.org / ). This density distribution can usually be regarded as a mixture of two similar normal distributions. Figure 2 ).

[0076] 5. Determine the genetic distance threshold

[0077] Considering that the density distribution obtained in step 4 exhibits a clear Gaussian mixture distribution, the Expectation-Maximization (EM) algorithm, implemented in the R package `mixtools`, is used to fit the density distribution obtained in step 4. This fits the mean and variance of the two normal distributions, and the cross boundary of the fitted normal distributions is selected as the hard threshold (i.e., the genetic distance threshold). In this invention, the hard threshold is specifically 10. -3 Mutations per nucleotide site (10) -3 per bp) as the dividing point between the two distributions ( Figure 3 ).

[0078] 6. Divide the initial ancestral haplogroups

[0079] The genetic distance matrix of each window output by PLINK v1.9 in step 4 is clustered using the average-linkage hierarchical clustering method to obtain the hierarchical clustering tree for each window.

[0080] Using the genetic distance hard threshold obtained in step 5 above as the height, the hierarchical clustering tree of each small window is further divided into several data groups to obtain the initial ancestral haplogroup.

[0081] 7. Determine the initial priority order of samples

[0082] For the initial ancestral haplotype obtained in 6 above, a greedy algorithm is used to obtain a sample subset that maximizes the diversity of ancestral haplotypes. In each round, the sample that shares the most ancestral haplotypes with other samples outside the subset is selected into the subset until most (by default, 95% or more) of the ancestral haplotypes are covered by the samples in the subset. Then, the samples not covered by the subset are added to the end of the subset priority in random order. The order in which the samples are added to the subset is the initial priority order of the samples.

[0083] 8. Determine the ancestral haplotype type of all windows

[0084] According to the haplotype of the sample to be tested, Zavitan is classified into the tetraploid wheat group, and Laolaihou, Hongchunmai, and Lovrin 10 are classified into the hexaploid wheat group. The semi-supervised mode is used to determine the ancestral haplotype type of all windows, and the priority order between the two groups is tetraploid wheat group > hexaploid wheat group.

[0085] Using the bubble-like algorithm, the initial ancestral haplotype obtained in 6 above and the initial priority order of the samples obtained in 7 above are used to obtain an ancestral haplotype file and a sample priority order file that record the ancestral haplotype type of all windows. After the above calculation, the priority order of the samples is changed so that the total number of changes in ancestral haplotype between adjacent windows is minimized.

[0086] Using the semi-supervised mode, the priority order of the samples is changed using the bubble-like algorithm so that the total number of changes in ancestral haplotype between adjacent windows is minimized. Starting from the sample with the highest priority, the priority of two adjacent samples is exchanged, and the ancestral haplotype is determined according to the priority order of the samples before and after the exchange, respectively. If the total number of changes is reduced, the exchange is retained, otherwise the priority order is restored. Until the exchange of sample order cannot reduce the total number of haplotype changes, the step is ended. Then, the priority order of the samples is combined according to the priority order between the groups to obtain a sample priority order file, and the ancestral haplotype type of all windows is determined according to this sample priority order to obtain an ancestral haplotype file.

[0087] After determining the priority order of the samples in each round, the method for determining the ancestral haplotype is as follows 1)-3):

[0088] (1) Assign the same new ancestral haplotype type to all windows of the sample with the highest priority;

[0089] (2) Starting from the sample with the second highest priority, for each window of each sample, if there is a sample with a higher priority that has the same initial ancestral haplotype as the sample in the window, the ancestral haplotype type of the sample in the window is determined to be the same as the ancestral haplotype type of the high-priority sample;

[0090] (3) From the second highest priority sample, for each small window of each sample, if there is no higher priority sample with the same initial ancestral haplotype as the sample in the small window, the ancestral haplotype type of the sample in the small window is determined as a new ancestral haplotype type.

[0091] After the ancestral haplotype type is determined, the total number of ancestral haplotype changes is calculated as follows: for each sample, if the ancestral haplotype types of adjacent small windows are different, the total number of ancestral haplotype changes is increased by 1.

[0092] 9. Adjust and optimize the ancestral source interval

[0093] The ancestral haplotype file obtained in the above 8 is corrected by using the Bayesian method to obtain the corrected ancestral haplotype file, which is the target ancestral haplotype file. The file contains the ancestral haplotype attribution state of each small window on the whole genome of the sample to be tested. The above correction fine-tunes the interval division and reduces the random errors caused by the hard threshold.

[0094] The above correction refers to the adjacent small window to correct the possible error distribution of haplotype under the hard threshold. For each small window, the possibility of error distribution is first estimated by the ratio between type I error (a) and type II error (β) of the current haplotype classification. The two types of errors are quantified by Gaussian distribution based on the genetic distance from the representative joiner. If the probability of type II error (β) is greater than a certain multiple (default is 10) of type I error (a), the haplotype of the small window is marked as correct and no change is made. Otherwise, consider the possibility of all haplotype types in the upstream and downstream 10 windows of this small window, calculate the product of the β value and the occurrence frequency of each haplotype type, and reassign the small window to the ancestral haplotype type with the maximum product value.

[0095] The corrected ancestral haplotype file can be visualized by using R language software (https: / / cran.r-project.org / ) to visualize the results of step 9 to obtain more intuitive distribution information of the same or different ancestral source interval on each chromosome (as shown in Figure 4 ).

[0096] Figure 4 As can be seen, the segments with the same color in different samples to be tested represent the same ancestral source, indicating that the samples have the same ancestral source in the segment, and the segment is named as the ancestral haplotype block.

[0097] Example 2, verification of the fine interval ancestral haplotype typing method of the wheat genome

[0098] For the segments with different ancestral haplotypes obtained in Example 1, the software designed genome-specific primers were used for amplification and amplification product sequencing verification.

[0099] 1. Primer design and synthesis

[0100] From the segments with different ancestral haplotypes obtained in Example 1, Figure 4 It can be seen that the wheat material Lovell No. 10 to be tested and Hongchunmai belong to different ancestral haplotypes in the region of chromosome 6A (chr6A) 400Mbp-440Mbp, and the wheat material Laolaihou and Hongchunmai belong to the same ancestral haplotype. Taking the region as an example, the same ancestral source interval is designed as follows:

[0101] F: TAGATCGATGCTCTTGAGGC (sequence 1); R: AGAAGTTCGAAGTGGCTCTC (sequence 2).

[0102] The amplification primers were designed by online software Primer 3.0, wherein the length of the primers was 18-24bp, the GC content was 40%-60%, the annealing temperature was 54-60℃, and the amplification product size of the amplification primers was 100-500bp. The primers were synthesized by Beijing Genesee Biotechnology Co., Ltd.

[0103] 2. PCR amplification and product sequencing

[0104] The CTAB method was used to extract the DNA of the wheat materials Laolaihou, Hongchunmai and Lovell No. 10 to be tested, wherein Laolaihou was used as a control for the method in Example 1.

[0105] The above genome-specific primers were used for fragment amplification, and the obtained PCR products were sent to Beijing Norsen Genomic Research Center Co., Ltd. for Sanger sequencing verification of the difference in nucleotide sequence in the selected DNA segment.

[0106] The above PCR reaction system was 20μL, including 10μL of 2X M5 HiPer plus Taq HiFi PCR mix, 1μL of forward and reverse primers (10μmol L –1 ) each, 150ng μL –1 template DNA 2μL, and ddH2O was added to 20μL. The PCR amplification program was 95℃ for 3min; 95℃ for 30s, 56-57.4℃ for 30-60s (depending on the primer annealing temperature and target sequence), 72℃ for 2min, 35 cycles; 72℃ for 5min.

[0107] The PCR product sequencing results are as follows: Figure 5As shown, no single nucleotide polymorphisms exist between the wheat materials Laihao and Hongchunmai in the selected interval of the same ancestral origin of the two varieties analyzed by the method in Example 1, while the control material Lovrin 10 has partial single nucleotide differences with the two varieties in the selected interval.

[0108] The above results show that the method in Example 1 can effectively screen or identify the blocks of the test varieties with the same ancestral haplotype, or the blocks of different ancestral haplotypes.

Claims

1. A method for determining the ancestral haplogroup to which different wheat sample genome blocks belong, comprising the following steps: 1) Obtain whole-genome resequencing data for each wheat sample to be tested; 2) Filter the whole genome resequencing data of each wheat sample to be tested to achieve quality control and obtain filtered data; then compare the filtered data with the reference wheat genome sequence to obtain the DNA sequence alignment file of each sample to be tested; filter the DNA sequence alignment file of each sample to be tested again to obtain the filtered DNA sequence alignment file of each sample to be tested. 3) Perform SNP and CNV variant identification on the DNA sequence alignment files of all the samples to be tested after filtering, and obtain the filtered VCF file containing all SNP site information and the BED file containing CNV variant interval information for each sample to be tested. 4) Take the filtered VCF file containing all SNP locus information and the BED file containing CNV variation interval information for each sample to be tested, find and count the number of differential SNPs in the small window where neither of the two samples is a CNV variation interval, divide the number of each SNP by the corresponding small window length to get the density of each SNP, and obtain the genetic distance matrix file that records the density of SNPs that are different between the two varieties in the small window; then plot the distribution of SNP polymorphism site density based on the SNP density of all samples. 5) After fitting the density distribution of the SNP polymorphic sites to obtain a normal distribution, determine the genetic distance hard threshold based on the boundary; the fitting algorithm is the expectation-maximization algorithm; The hard threshold is determined using the cross boundary of the fitted normal distribution as the hard threshold; 6) Cluster the genetic distance matrix of each small window obtained in step 4) to obtain a hierarchical clustering tree sample for each small window; then use the genetic distance hard threshold obtained in step 5) as the height to divide the hierarchical clustering tree sample of each small window to obtain the initial ancestral haplogroup. 7) Calculate the initial ancestral haplogroups to obtain a sample subset that maximizes the diversity of ancestral haplogroups, and then add the samples that were not selected into the subset in a random order; use the order in which the samples are selected into the subset as the initial priority order of the samples. 8) Calculate the initial ancestral haplogroup obtained in step 6) and the sample priority order file obtained in step 7) to obtain the sample priority order file and the ancestral haplogroup file; the purpose of the calculation is to minimize the total number of ancestral haplogroup changes between adjacent windows.

2. The method according to claim 1, characterized in that: The method further includes the following step 9): correcting the ancestral haplogroup file to obtain the corrected ancestral haplogroup file.

3. The method according to claim 2, characterized in that: The method further includes the following step 10): the corrected ancestral haplogroup file is visualized to obtain a visualization map of the ancestral haplogroup attribution status of each segment on the whole genome of each test sample; based on the segments with the same color in different test samples in the visualization map, the segment is determined to be the ancestral haplogroup block.

4. The method according to claim 1 or 2, characterized in that: In step 1), the wheat sample is a variety or strain, or a closely related species or a wild species.

5. The method according to claim 1 or 2, characterized in that: In step 4), the length of the small window in the CNV variant identification is 1-5 Mbp.

6. The method according to claim 1 or 2, characterized in that: In step 8), the calculation method is any one of the following two modes: 8-1) Unsupervised mode, where the priority order of samples is completely determined by a bubble-like algorithm; 8-2) Semi-supervised mode, where the priority order of samples is determined by a bubble-like algorithm within a predefined sample group.

7. The method according to claim 2, characterized in that: In step 9), the correction is performed using a Bayesian method.

Images