A method and application of InDel background molecular marker design based on resequencing

Abstract

The present invention discloses a method and application of designing InDel background molecular markers based on resequencing, including a data preprocessing module, a variant site screening module, a primer design module and a background molecular marker development module. The present invention realizes the identification and screening of InDel markers in the whole genome based on resequencing data, and performs primer design and background molecular marker development, and has the advantages of flexible operation, high degree of automation and high efficiency in batches. The present invention designs background molecular markers based on the resequencing data of the sample. Compared with the traditional design of primers based on the reference genome, the primers have stronger specificity, higher density and more quantity. The present invention innovatively uses windowing to screen InDel mutations, unifies that only one InDel marker is allowed to exist in each window, reduces the uncertainty of directly finding variant sites based on the reference genome, and ensures the reliability and genotyping ability of background molecular markers.

Description

Technical Field

[0001] The present invention relates to the technical field of sequencing, and in particular to a method and application of designing InDel background molecular markers based on resequencing. Background Art

[0002] The collection, identification and evaluation of crop germplasm resources are important measures to maintain species genetic diversity and ensure ecological balance. DNA insertion / deletion (InDel) markers have the characteristics of simple banding, high density, easy detection and accurate results, and are very suitable as background molecular markers for genetic identification and evaluation of germplasm resources. Next generation sequencing (NGS), also known as resequencing, can identify genetic variation within the genome, providing valuable information for breeding and crop genetic improvement. In summary, the development of corresponding background molecular markers for genomic genetic variation is conducive to the rapid identification and evaluation of germplasm resources at the gene level, improving work efficiency and reducing costs.

[0003] Traditional marker development is based on primer design based on a reference genome sequence. However, due to the lack of individual genetic variation information, the development of primers based on traditional methods is very blind, usually requires experimental verification, and has low usability. In addition, the development of markers based on genetic variation information requires a lot of manual screening and repeated confirmation in the early stage, which greatly increases the time and money costs. Summary of the invention

[0004] The object of the present invention is to provide a method and application of designing InDel background molecular markers based on resequencing to solve the problems raised in the above-mentioned background technology.

[0005] To achieve the above objectives, the present invention provides the following technical solutions: a device for designing InDel background molecular markers based on resequencing, comprising a data preprocessing module, a variant site screening module, a primer design module and a background molecular marker development module.

[0006] Preferably, a method for designing a device for InDel background molecular markers based on resequencing comprises the following steps:

[0007] A. Sequencing data preprocessing:

[0008] A1. The fastq file of the raw sequencing data is named $samplename.R1 / 2.fastq.gz, where samplename is the sequencing sample name. Use ls and awk to write the sample name into a text file, with each line representing a sequencing sample.

[0009] A2. Use quality control software to perform quality detection on the original sequencing data, and remove adapters and low-quality reads;

[0010] A3. Use alignment software to align the re-sequencing data after quality control to the reference genome to obtain a bam format file containing information such as alignment positions, alignment quality, and linear relationship information with the reference genome;

[0011] A4. Use a repeat sequence detection tool to delete PCR repeats in the alignment results;

[0012] A5. Set the sample sequencing depth threshold and alignment rate threshold, and delete samples with sequencing depth or alignment rate below the threshold;

[0013] A6. Use genomic analysis tools to identify genomic genetic variations to obtain a set of genetic variations for each sample. If there are multiple genomic re-sequencing samples, repeat the above steps and then merge them;

[0014] A7. Extract the InDel dataset from the genomic genetic variation dataset in A6, save it in a VCF format file, and perform filtering;

[0015] A8. Add the ID of the InDel in the format of "ChrN__POS", such as "Chr5__21752";

[0016] B. InDel variant site screening:

[0017] B1. Convert the VCF file into a 012 genotype matrix, that is, use "0", "1", and "2" to represent the three genotypes of "0 / 0", "0 / 1", and "1 / 1" respectively. The InDel012 matrix includes information such as chromosome, start position, InDel number, reference sequence, variant sequence, reference sequence length, and variant sequence length.

[0018] B2. Use windowing to screen InDel markers;

[0019] B3. Set the minimum fragment size min_size and maximum fragment size max_size of the InDel variation, screen the variant sequence length altL, and select markers where min_size < altL < max_size;

[0020] C. InDel primer development:

[0021] C1. Extract the upstream and downstream sequences of the InDel as templates for primer design;

[0022] C2. Generate a primer design input file in Boulder-IO format, including InDel number, primer design parameters, template sequence, target sequence position, PCR amplification fragment size, and SNP variation coordinates that need to be avoided. After the template sequence uses regular expressions to extract variation coordinates, the annotations need to be deleted and then input into the Boulder-IO file;

[0023] C3. Use Pimer3 software to batch design InDel primers;

[0024] D. Background Molecular Marker Development:

[0025] D1. Extract the sequence information of InDel primers and generate a fasta format file for each pair of InDel primers;

[0026] D2. Use MFEprimer software to detect the PCR amplification specificity of InDel primers;

[0027] D3. Extract primer amplification sequence information, retain InDel primers that can specifically amplify the target fragment, and obtain candidate InDel background markers;

[0028] D4. Count the number of InDels in the amplification interval of the primers obtained in D3, and retain only the primers with only one InDel in the amplification interval in all sequencing materials to obtain the InDel background molecular marker;

[0029] D5. Using the genotype data, perform K-means clustering on the candidate background molecular markers obtained in D4 based on the specified number of clusters;

[0030] D6. Use R software to randomly sample the InDels of each cluster to obtain the core InDel background molecular markers.

[0031] Preferably, the filtering method in step A7 is as follows:

[0032] 1) Set a threshold for the number of supported reads and delete InDels that are not within the threshold range;

[0033] 2) Only retain biallelic InDels;

[0034] 3) Set the missing threshold of InDel genotype and delete the InDels whose genotype missing exceeds the threshold in all samples;

[0035] 4) Set the InDel minor allele frequency threshold and delete the InDels whose genotype proportion is less than the threshold.

[0036] Preferably, the marking method in step B2 is as follows:

[0037] 1) Set the window size and step size;

[0038] 2) Count the InDels in each window and keep the windows with only one InDel;

[0039] 3) Extract the InDel that meets the windowing condition and generate a new InDel012 matrix.

[0040] Preferably, the setting method of step C1 is:

[0041] 1) Specify the length ws of the upstream and downstream sequences to be extracted, with the optimal parameter being ws = 150 bp;

[0042] 2) The starting and ending positions of the primer template for fragment insertion on the genome are: start = POS–ws, end = POS+ws; the starting and ending positions of the primer template for fragment deletion on the genome are: start = POS–ws, end = POS+ws+altL;

[0043] 3) The target sequence for the fragment insertion primer amplification is: start = end = ws + 1; the target sequence for the fragment deletion primer amplification is: start = ws + 1, end = ws + altL + 1;

[0044] 4) Mark the variation type in the template sequence, use square brackets to mark InDel and parentheses to mark SNP.

[0045] Preferably, the sequencing raw data quality control software used in step A2 is FASTP, the sequence comparison software used in step A3 is BWA, the PCR duplicate deletion software used in step A4 is PICARD, the sequencing depth and alignment rate statistics software used in step A5 is SAMTOOLS, the genome analysis tool used in step A6 is GATK, and the InDel quality control software used in step A7 is GATK and VCFTOOLS; and step B1 uses VCFTOOLS and BCFTOOLS to generate the InDel012 matrix.

[0046] Compared with the prior art, the present invention has the following beneficial effects:

[0047] (1) The present invention realizes the identification and screening of InDel markers on the whole genome based on resequencing data, and performs primer design and background molecular marker development, which has the advantages of flexible operation, high degree of automation and high efficiency in batches.

[0048] (2) The present invention designs background molecular markers based on the resequencing data of the sample. Compared with the traditional primer design based on the reference genome, the primers have stronger specificity, higher density and greater quantity.

[0049] (3) The present invention innovatively uses windows to screen InDel variants, unifying that only one InDel marker is allowed to exist in each window, reducing the uncertainty of directly searching for variant sites based on the reference genome, and ensuring the reliability of background molecular markers and genotyping capabilities.

[0050] (4) The present invention integrates data quality control, InDel screening, primer design, and background molecular marker development steps to achieve automated operations from sequencing data to primers and then to background molecular marker sets. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 It is a schematic diagram of the main process of the present invention;

[0052] Figure 2 This is a schematic diagram of window screening of InDel markers of the present invention;

[0053] Figure 3 Schematic diagram for background molecular marker selection;

[0054] Figure 4 This is the core background molecular marker genotype map obtained by the present invention. DETAILED DESCRIPTION

[0055] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0056] It should be noted that the structures, proportions, sizes, etc. illustrated in the drawings in this specification are only used to match the contents disclosed in the specification so that people familiar with this technology can understand and read them. They are not used to limit the conditions under which this application can be implemented, so they have no substantive technical significance. Any structural modification, change in proportion or adjustment of size should still fall within the scope of the technical content disclosed in this application without affecting the effects and purposes that can be achieved by this application.

[0057] See also Figure 1-4 The present invention provides a technical solution: an InDel background molecular marker design device based on resequencing, comprising a data preprocessing module, a variant site screening module, a primer design module and a background molecular marker development module.

[0058] In the present invention, a method for designing an InDel background molecular marker device based on resequencing is provided, and the design method comprises the following steps:

[0059] A. Preprocessing of sequencing data:

[0060] A1. The naming of the original sequencing data fastq file is $samplename.R1 / 2.fastq.gz, where samplename is the name of the sequencing sample. Use ls and awk to write the sample names into a text file, with each line representing a sequencing sample.

[0061] A2. Use a quality control software to perform quality detection on the original sequencing data, removing adapters and low-quality reads.

[0062] A3. Use an alignment software to align the quality-controlled re-sequencing data to the reference genome, obtaining a bam format file containing information such as alignment positions, alignment quality, and the linear relationship with the reference genome.

[0063] A4. Use a repetitive sequence detection tool to delete PCR duplicates in the alignment results.

[0064] A5. Set the sample sequencing depth threshold and alignment rate threshold, and delete samples with a sequencing depth or alignment rate below the threshold.

[0065] A6. Use genomic analysis tools to identify genomic genetic variations, obtaining a set of genetic variations for each sample. If there are multiple genomic re-sequencing samples, repeat the above steps and then merge them.

[0066] A7. Extract the InDel dataset from the genomic genetic variation dataset in A6, save it in a VCF format file, and perform filtering.

[0067] A8. Add the ID of the InDel in the format of "ChrN__POS", such as "Chr5__21752".

[0068] B. Screening of InDel variant sites:

[0069] B1. Convert the VCF file into a 012 genotype matrix, using "0", "1", and "2" to represent the three genotypes of "0 / 0", "0 / 1", and "1 / 1" respectively. The InDel012 matrix includes information such as chromosome, start position, InDel number, reference sequence, variant sequence, reference sequence length, and variant sequence length.

[0070] B2. Use windowing to screen for InDel markers.

[0071] B3. Set the minimum fragment size min_size and maximum fragment size max_size of the InDel variation, screen the variant sequence length altL, and select markers where min_size < altL < max_size.

[0072] C. InDel primer development:

[0073] C1. Extract the upstream and downstream sequences of InDel and use them as templates for primer design;

[0074] C2. Generate a primer design input file in Boulder-IO format, including InDel number, primer design parameters, template sequence, target sequence position, PCR amplification fragment size, and SNP variation coordinates that need to be avoided. After the template sequence uses regular expressions to extract variation coordinates, the annotations need to be deleted and then input into the Boulder-IO file;

[0075] C3. Use Pimer3 software to batch design InDel primers;

[0076] D. Background Molecular Marker Development:

[0077] D1. Extract the sequence information of InDel primers and generate a fasta format file for each pair of InDel primers;

[0078] D2. Use MFEprimer to detect the PCR amplification specificity of InDel primers;

[0079] D3. Extract primer amplification sequence information, retain InDel primers that can specifically amplify the target fragment, and obtain candidate InDel background markers;

[0080] D4. Count the number of InDels in the amplification interval of the primers obtained in D3, and retain only the primers with only one InDel in the amplification interval in all sequencing materials to obtain the InDel background molecular marker;

[0081] D5. Using the genotype data, perform K-means clustering on the candidate background molecular markers obtained in D4 based on the specified number of clusters;

[0082] D6. Use R software to randomly sample the InDels of each cluster to obtain the core InDel background molecular markers.

[0083] In the present invention, the filtering method in step A7 is as follows:

[0084] 1) Set a threshold for the number of supported reads and delete InDels that are not within the threshold range;

[0085] 2) Only retain biallelic InDels;

[0086] 3) Set the missing threshold of InDel genotype and delete the InDels whose genotype missing exceeds the threshold in all samples;

[0087] 4) Set the InDel minor allele frequency threshold and delete the InDels whose genotype proportion is less than the threshold.

[0088] In the present invention, the marking method of step B2 is as follows:

[0089] 1) Set the window and step size;

[0090] 2) Count the InDels in each window and keep the windows with only one InDel;

[0091] 3) Extract the InDel that meets the windowing condition and generate a new InDel012 matrix.

[0092] In the present invention, the setting method of step C1 is:

[0093] 1) Specify the length ws of the upstream and downstream sequences to be extracted, with the optimal parameter being ws = 150 bp;

[0094] 2) The starting and ending positions of the primer template for fragment insertion on the reference genome are: start = POS–ws, end = POS+ws; the starting and ending positions of the primer template for fragment deletion on the genome are: start = POS–ws, end = POS+ws+altL, POS is the starting position of the variation on the chromosome, and altL is the length of the deletion;

[0095] 3) The target sequence for the fragment insertion primer amplification is: start = end = ws + 1; the target sequence for the fragment deletion primer amplification is: start = ws + 1, end = ws + altL + 1;

[0096] 4) Mark the variation type in the template sequence, use square brackets to mark InDel and parentheses to mark SNP.

[0097] In the present invention, step A2 uses the sequencing raw data quality control software FASTP, step A3 uses the sequence comparison software BWA, step A4 PCR duplicate deletion software PICARD, step A5 sequencing depth and alignment rate statistics software SAMTOOLS, step A6 genome analysis tool GATK, step A7 InDel quality control software GATK and VCFTOOLS; step B1 uses VCFTOOLS and BCFTOOLS to generate the InDel012 matrix, FASTP, BWA, PICARD, SAMTOOLS, GATK, VCFTOOLS and BCFTOOLS are software and have no Chinese names, and are directly expressed in English in the industry.

[0098] The sequencing depth threshold and alignment rate threshold of A5 of the present invention are designed by the user. For example, if they are set to 5 and 70, it means that the sequencing samples with an average sequencing depth of less than 5 or an alignment rate of less than 70% are deleted; the InDel support number threshold of A7 of the present invention is designed by the user, such as 5-40, which means that only InDels with reads support numbers greater than 5 and less than 40 are retained; the InDel genotype deletion threshold of A7 of the present invention is designed by the user, such as 0.9, which means that the InDel genotype deletion rate in all materials exceeds 10%, and the InDel is deleted; the minor allele frequency threshold of A7 of the present invention is designed by the user, such as 0.05, which means that the InDel with a minor allele frequency of less than 0.05% is deleted; the window and step size of B2 of the present invention are designed by the user, such as 300 and 50, and the InDels in the window are counted according to the window of windows=300bp and step=50bp; B3 of the present invention The size of the InDel variant fragment is designed by the user, such as 20-100, which means that the InDel with a fragment of 20 to 100 bp is screened; the interval of the upstream and downstream sequences of InDel extracted by C1 of the present invention is designed by the user, such as 150, which means that the 150 bp sequence upstream and downstream of InDel is extracted as the template sequence for primer design; the primer design parameters in C2 of the present invention are designed by the user, such as 182022100-200, which means that the primer length is between 18-22 bp, the optimal primer length is 20 bp, and the amplified fragment size is 100-200 bp; the number of clusters in D5 of the present invention is designed by the user, such as 100, which means that the number of K-mean clusters is 100. All the above parameters that can be designed by the user have default parameters, and the user only needs to specify the corresponding parameters according to the needs.

[0099] Embodiment 1:

[0100] Resequencing data preprocessing:

[0101] Based on 157 core germplasms of South China indica rice whose genomes have been resequenced, the present invention is used to design core background molecular markers within the genome.

[0102] The sequencing data file is named $samplename.R1 / 2.fastq.gz, where samplename is the sequencing sample name. $samplename is written into the samplename.txt text file, and the self-written program 1fastp_batch.sh is run to perform quality control on the original resequencing data. The quality control software is FASTP, and the parameters are set to: -q 20-u 40-l 50-n 10-f 12-F 12. Since the sequencing data of each sample are independent of each other, the samplename.txt can be split into multiple files to achieve a similar effect of multi-threaded parallel computing.

[0103] Run the self-written program 2bwa_men_map.sh to align the resequencing data after quality control to the reference genome and obtain a sam file containing the alignment results. Use SAMTOOLS software to sort the sam file and convert it into a bam file with smaller storage space. The alignment software used is BWA, and the reference genome is Nipponbare MSU7.0.

[0104] PCR duplicates were removed from the bam file using the MarkDuplicates module of PICARD with the parameter set to REMOVE_DUPLICATES = true and other parameters as default values.

[0105] The sequencing depth and alignment rate were calculated using SAMTOOLS and awk. The sequencing depth of 157 South China indica rice samples ranged from 2.6 to 42.5×, with an average sequencing depth of 30.7× and an alignment rate greater than 70%. After deleting two samples with a sequencing depth lower than 5×, the remaining sequencing samples were subjected to subsequent analysis.

[0106] GATK software was used to extract SNP and Indel site information from bam files, and GVCF files were generated for each sample. GATK software was used to merge all sample GVCF files. The genetic variation information was extracted from the merged GVCF files using the GATK GenotypeGVCFs parameter, and a total of 20,740,579 SNPs and 5,726,524 InDels were identified. GATK and VCFTOOLS were used to filter the identified InDels and SNPs. GATK parameters were set as follows: --filter-expression "QD<2.0||FS>60.0||MQ<40.0||SOR>3.0||MQRankSum<-12.5||ReadPosRankSum<-8.0" --genotype-filter-expression "DP<5||DP>40" --restrict-alleles-to BIALLELIC; VCFtools parameters were set as follows: --max-missing 0.9 --maf 0.05. After filtering, VCF format files containing 2,246,018 and 542,957 high-quality SNPs and InDels were obtained, respectively. Run the self-written program 3VCF_add_id.pl to add ID information to the VCF file. The ID format is "ChrN__POS", such as "Chr5__21752".

[0107] Example 2: InDel marker screening

[0108] VCFTOOLS was used to convert the InDel VCF file into a matrix with "0", "1" and "2" representing the three genotypes of "0 / 0", "0 / 1" and "1 / 1", respectively, to obtain a .012 file containing the genotype, a .indv file containing the sample name, and a .pos file containing the InDel position. BCFTOOLS (parameters: query-f'%CHROM\t%POS\t%ID\t%REF\t%ALT\n') was used to generate an ALT.bed file containing the chromosome position, reference sequence and variant sequence.

[0109] The self-written program 4vcf012matrix.R was used to merge the information, generate the InDel012 matrix, and screen the markers. Given the parameters: 1) output file prefix; 2) chromosome interval, such as "Chr5:26545-589242", if the parameter is "all", all InDels will be generated into a 012 matrix; 3) window size; 4) step size; 5) minimum variant fragment length; 6) maximum variant fragment length, such as Rscript vcf012ma trix.R DEC31 all 3001020100). This self-written program implements the following functions:

[0110] 1. Merge the .012, .indv, .pos and ALT.bed files to generate the InDel012 matrix. If a chromosome interval is specified, the InDel012 matrix within the specified interval is generated. Write the InDel012 matrix of this step into a TXT file named $prefix.genotype.raw.txt, where prefix is ​​the user-defined output file prefix.

[0111] 2. If Figure 2 As shown, the InDel012 matrix is ​​windowed using the window (window) and step (step), and the number of InDels in each window is calculated. The window with only one InDel is selected to generate a new InDel012 matrix; in Example 2, window = 300bp, step = 10bp. After windowing and screening, a total of 256,387 InDels are retained for subsequent analysis;

[0112] 3. InDel fragment size screening. In Example 2, the InDel fragment size was selected to be 20-100 bp, and a total of 96,387 InDels were obtained. The InDel012 matrix of this step was written into a TXT file named $prefix.genotype.txt, where prefix is ​​the user-defined output file prefix.

[0113] Example 3: InDel primer design

[0114] Run the self-written program 5InDel_Primer.py to design InDel primers, and give the following parameters: 1) -g: specify the reference genome; 2) -p: specify the output file prefix; 3) -s: enter a value to specify how many bp sequences upstream and downstream of InDel are extracted for the primer design template sequence, and the default value is 150; 4) -opt: enter a value to specify the optimal primer length, and the default value is 20; 5) -min: enter a value to specify the minimum primer length, and the default value is 18; 6) -max: enter a value to specify the maximum primer length, and the default value is 22; -ps: enter a value to specify the size range of the primer amplification fragment, and the default value is 100-200. Among the parameters, -g and -p are required parameters, and other parameters can call the default values. The program calls the $prefix.genotype.txt file generated in Example 2 by default. The self-written program realizes the following functions:

[0115] 1. Automatically determine whether there is an index file for the reference genome. If not, call MFEprimer software to generate a reference genome index file;

[0116] 2. Obtain the target sequence coordinates for primer design. Based on $prefix.genotype.txt, write the position, ID, primer template sequence coordinates and primer target sequence coordinates of the InD el marker on the chromosome into the forprimer_withposID.txt file. In the embodiment, the target sequence for fragment insertion primer amplification is: start=end=151; the target sequence for fragment deletion primer amplification is: start=151, end=151+altL, where altL is the length of the deletion fragment. Write the position, ID and coordinates of the SNP marker on the chromosome into the SNP.bed file;

[0117] 3. Use the POS information in the $prefix.genotype.txt file and the value specified by -s to generate the subseq.bed file containing the chromosome position information of the primer design template sequence. In the embodiment, specify -s=150, then the starting and ending positions of the primer template for fragment insertion on the genome are: start=POS-150, end=POS+150; the starting and ending positions of the primer template for fragment deletion on the genome are: start=POS-150, end=POS+150+altL, POS is the starting position of the variation on the chromosome, and altL is the length of the deletion fragment;

[0118] 4. Based on the subseq.bed file, call SEQTK software to extract the template sequence for primer design from the reference genome and write it to the file tmp.fa;

[0119] 5. Call the fx2tab module of SEQKIT software to convert tmp.fa into a tmp2.fa file in tabular form;

[0120] 6. Annotate the variant sites in the primer design template sequence of the tmp2.fa file, where InDel sequences are marked with "[]" and SNP variants are marked with "()";

[0121] 7. Merge tmp2.fa and forprimer_withposID.txt to generate withseq.txt file;

[0122] 8. Generate the standard format file primer3_Boulder-IO_input.txt required for Pimer3 software input based on withseq.txt. First extract the mutation coordinates of the primer design template sequence, and then input the sequence with the mutation site annotation deleted into the primer3_Boulder-IO_input.txt file;

[0123] 9. Call Primer3 software to generate InDel primers in batches, organize the software output results into a table, and generate the primer_0_info.txt file.

[0124] Example 4: Development of InDel background molecular markers

[0125] The operation process of this embodiment is as follows Figure 3 As shown, the first line runs the self-written program 6PrimerAmpDetecti on.py to perform InDel primer specific detection, with given parameters: 1) -g: specify reference genome; 2) -p: specify output file prefix. The self-written program runs to achieve the following functions:

[0126] 1. Based on primer_0_info.txt, extract the InDel primer sequence and generate a / temp / $ID.fa file for each pair of primers. The ID is the ID of the InDel marker.

[0127] 2. Call the spec module of MFEprimer software to detect the specificity of the primers in 1 on the reference genome;

[0128] 3. Use regular expressions to organize the output results of MFEprimer, screen primers with only one potential amplification fragment, and write the InDel ID and potential amplification range into the $prefix_specific_primer.txt file.

[0129] Run the self-written program 7primer_info_combine.R to integrate the primer_0_info.txt, $prefix_specific_primer.txt and $prefix.genotype.txt files to obtain the InDel primer set for specific amplification and generate the $prefix.vcfprimer.txt file, as shown in the following table:

[0130]

[0131] Run the self-written program 8kmeans_slection.R to perform cluster analysis on the InDe l markers in $prefix.vcfprimer.txt, given the parameters 1) output file prefix; 2) the expected number of kmean clusters. The self-written program implements the following functions:

[0132] 1. Based on the $prefix.vcfprimer.txt and $prefix.genotype.raw.txt files, call BEDTOOLS to analyze the number of InDels in the primer amplification interval, retain the primers with only one InDel in the amplification interval, and obtain all qualified InDel background markers;

[0133] 2. Based on the genotype information in $prefix.vcfprimer.txt, call R's kmeans function to perform K-means cluster analysis;

[0134] 3. Call the R package sampling and use stratified sampling to select an InDel primer from each kmeans cluster to form the InDel core background molecular marker.

[0135] Based on the genotype data of 155 core germplasms of indica rice in South China, 3161 InDel markers were clustered by K-mean, and 100 pairs of core background molecular markers were obtained. Figure 4 As shown, the use of 100 pairs of core markers can be used to distinguish the core germplasm of indica rice in South China, as well as evaluate genetic diversity and kinship.

[0136] In summary, the present invention realizes the identification and screening of InDel markers in the whole genome based on resequencing data, and performs primer design and background molecular marker development, which has the advantages of flexible operation, high degree of automation and high efficiency in batches; the present invention designs background molecular markers based on the resequencing data of the sample, and compared with the traditional design of primers based on the reference genome, the primers have stronger specificity, higher density and more quantity; the present invention innovatively uses windowing to screen InDel mutations, unifies that only one InDel marker is allowed in each window, reduces the uncertainty of directly searching for mutation sites based on the reference genome, and ensures the reliability and genotyping ability of background molecular markers; the present invention integrates data quality control, InDel screening, primer design, and background molecular marker development steps, and realizes automated operations from sequencing data to primers and then to background molecular marker sets.

[0137] It should also be noted that, in this specification, relational terms such as first and second, etc. are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "comprise", "include" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, an element defined by the statement "comprises a ..." does not exclude the presence of other identical elements in the process, method, article or device including the element.

Claims

1. A device for designing InDel background molecular markers based on resequencing, characterized in that: It includes a data preprocessing module, a variant site screening module, a primer design module, and a background molecular marker development module; The design method includes the following steps: A. Sequencing data preprocessing: A1. The naming of the original sequencing data fastq file is $samplename.R1 / 2.fastq.gz, where samplename is the name of the sequencing sample. Use ls and awk to write the sample names into a text file, with each line representing a sequencing sample; A2. Use quality control software to perform quality detection on the original sequencing data, removing adapters and low-quality reads; A3. Use alignment software to align the quality-controlled resequencing data to the reference genome to obtain a bam format file containing alignment positions, alignment quality, and information on the linear relationship with the reference genome; A4. Use a repeat sequence detection tool to delete PCR repeats in the alignment results; A5. Set the sample sequencing depth threshold and alignment rate threshold, and delete samples with a sequencing depth or alignment rate below the threshold; A6. Use genomic analysis tools to identify genomic genetic variations, obtain the genetic variation set of each sample. If there are multiple genomic resequencing samples, repeat the above steps and then merge them; A7. Extract the InDel dataset from the genomic genetic variation dataset in A6, save it in a VCF format file, and perform filtering; A8. Add the ID of the InDel in the format of "ChrN__POS”; B. InDel variant site screening: B1. Convert the VCF file into a 012 genotype matrix, using "0", "1", and "2" to represent the three genotypes of "0 / 0", "0 / 1", and "1 / 1" respectively. The InDel012 matrix includes chromosome, start position, InDel number, reference sequence, variant sequence, reference sequence length, and variant sequence length information; B2. Use windowing to screen for InDel markers; B3. Set the minimum fragment size min_size and maximum fragment size max_size of the InDel variant, screen the variant sequence length altL, and select markers where min_size < altL < max_size; C. InDel primer development: C1. Extract the upstream and downstream sequences of the InDel for use as templates for primer design; C2. Generate a primer design input file in Boulder-IO format, including InDel number, primer design parameters, template sequence, target sequence position, PCR amplification fragment size, and SNP variant coordinates that need to be forced to avoid. After extracting the variant coordinates from the template sequence using regular expressions, the annotations need to be deleted before inputting into the Boulder-IO file; C3. Use Pimer3 software to perform batch design of InDel primers; D. Background molecular marker development: D1. Extract the sequence information of the InDel primers and generate a fasta format file for each pair of InDel primers; D2. Use MFEprimer software to detect the PCR amplification specificity of the InDel primers; D3. Extracting primer amplification sequence information, retaining InDel primers that can specifically amplify the target fragment, and obtaining candidate InDel background markers; D4. Count the number of InDels in the amplification interval of the primers obtained in D3, and retain only the primers with only one InDel in the amplification interval in all sequencing materials to obtain the InDel background molecular marker; D5. Using the genotype data, perform K-means clustering on the candidate background molecular markers obtained in D4 based on the specified number of clusters; D6. Using R software, randomly sample the InDels of each cluster to obtain the core InDel background molecular markers; The filtering method in step A7 is as follows: 1) Set a threshold for the number of supported reads and delete InDels that are not within the threshold range; 2) Only retain biallelic InDels; 3) Set the missing threshold of InDel genotype and delete the InDels whose genotype missing exceeds the threshold in all samples; 4) Set the InDel minor allele frequency threshold and delete the InDels whose genotype proportion is less than the threshold; The marking method of step B2 is as follows: 1) Set the window size and step size; 2) Count the InDels in each window and keep the windows with only one InDel; 3) Extract the InDel that meets the windowing conditions and generate a new InDel012 matrix; The setting method of step C1 is: 1) Specify the length ws of the upstream and downstream sequences to be extracted; 2) The starting and ending positions of the primer template for fragment insertion on the genome are: start = POS – ws, end = POS + ws; the starting and ending positions of the primer template for fragment deletion on the genome are: start = POS – ws, end = POS + ws + altL; 3) The target sequence for the fragment insertion primer amplification is: start = end = ws + 1; the target sequence for the fragment deletion primer amplification is: start = ws + 1, end = ws + altL + 1; 4) Mark the variant type in the template sequence, use square brackets to mark InDel and parentheses to mark SNP; The sequencing raw data quality control software used in step A2 is FASTP, the sequence comparison software used in step A3 is BWA, the PCR duplicate removal software used in step A4 is PICARD, the sequencing depth and alignment rate statistics software used in step A5 is SAMTOOLS, the genome analysis tool used in step A6 is GATK, and the InDel quality control software used in step A7 is GATK and VCFTOOLS; step B1 uses VCFTOOLS and BCFTOOLS to generate the InDel012 matrix.

Images