Mutation detection method based on end-to-end assembly genome

Abstract

The invention discloses a variation detection method based on an end-to-end assembled genome, and relates to a variation detection method. The objective of the invention is to solve the problems of incomplete detection, inaccurate detection, incapability of identifying complex variation types and low genotyping accuracy existing in a current method for detecting structural variation based on splicing. The variation detection method based on the end-to-end assembly genome comprises the following steps: step 1, processing a variation signal to obtain a processed variation signal; 2, performing mixed sorting on the processed variation signals obtained in the step 1, performing spatial distribution clustering on the sorted variation signals, judging whether the variation signals subjected to spatial distribution clustering are potential complex variation or not, if yes, judging the type of the complex variation, outputting the complex variation, obtaining complex variation signals, and executing the step 3; if not, indicating that the signal is a simple variation signal, and executing the step 3; and step 3, carrying out genetic typing on the variation signal obtained in the step 2. The method is applied to the field of variation detection.

Description

Technical Field

[0001] The present invention relates to a method for detecting variation. Background Art

[0002] At present, the means of analyzing the gene sequence of an organism is to obtain DNA sequence fragments with high credibility through sequencing and quality control, and then analyze the obtained fragments according to the corresponding tasks through various tools and methods. The most important analysis is to obtain the differences from the known reference sequence, that is, the variation in the biological sense; the variation with longer length and more complex form is called structural variation, including insertion (INS), deletion (DEL), duplication (DUP), translocation (TRA) and inversion (INV) and their fusion forms. Existing evidence shows that structural variation has a huge impact on the characteristics of organisms, so it is very important to detect structural variation. At present, most of the common processes for detecting structural variation in the human genome are based on second-generation or third-generation sequencing data, but due to the improvement of splicing tools and technological advances, it is possible to generate end-to-end spliced ​​chromosomes for individual typing. Existing studies have also shown that structural variation detection based on spliced ​​genomes can discover variations that may not be found in the second and third generations and has advantages in detecting complex structural variations. However, the current methods for detecting structural variation based on splicing have problems such as incomplete detection, inaccurate detection, inability to identify complex variation types, and the accuracy of genotyping still needs to be improved.

[0003] Therefore, this method provides a set of processes and methods for detecting structural variations based on end-to-end splicing, which improves the accuracy, sensitivity and detection efficiency of existing genome structural variation detection. Summary of the invention

[0004] The purpose of the present invention is to solve the problems of incomplete detection, inaccurate detection, inability to identify complex variation types and low genotyping accuracy in the current methods based on splicing to detect structural variation, and to propose a variation detection method based on end-to-end genome assembly.

[0005] The specific process of the variation detection method based on end-to-end assembled genome is as follows:

[0006] Step 1: Process the variant signal to obtain the processed variant signal; the specific process is:

[0007] Step 11, obtain BAM format comparison information;

[0008] Step 12: Filter the BAM format comparison information obtained in step 11 to obtain filtered BAM format comparison information;

[0009] Step 13: extract the variation signal from the filtered BAM format comparison information;

[0010] Step 14: Cluster the extracted variant signals and output the clustering results;

[0011] Step 15: Filter all the output cluster results to obtain the final cluster result, that is, integrate the signals contained in each cluster into a structural variation signal for output;

[0012] Step 2: perform mixed sorting on the processed variation signals obtained in step 1, perform spatial distribution clustering on the sorted variation signals, and determine whether the variation signals after spatial distribution clustering are potential complex variations. If so, it is determined to be a complex variation signal, and a complex variation signal is obtained, and step 3 is executed; if not, it is determined to be a simple variation signal, and step 3 is executed. A simple variation signal is a variation signal of insertion, deletion, inversion, translocation, and duplication;

[0013] Step 3: Perform genotyping on the variant signal obtained in step 2.

[0014] The beneficial effects of the present invention are:

[0015] The present invention provides a set of processes and methods for a variation detection method based on end-to-end genome assembly, which divides the variation detection process into three parts: a variation signal processing part, a complex structural variation detection part and a genotyping part, so as to realize efficient and reliable detection of structural variation based on splicing genome. In variation signal processing, a comparison information filtering method that maximizes non-overlapping comparison information is screened out to reduce false positives in variation detection; in variation signal extraction, a rule scheme for detecting basic structural variation is given, and different threshold parameters are set for different variation types, so as to extract variation signals more accurately and comprehensively; in variation signal clustering and filtering, an unsupervised clustering and filtering method based on feature similarity is proposed, which measures the similarity of variation signals through a similarity function to achieve clustering of similar signals in each variation type, forming an average variation to achieve a filtering effect; at the same time, abnormal variation aggregation signals are discarded to improve the credibility of the detected variation signals. In complex structural variation detection, by identifying potential complex variation sites and then performing complex structural variation fusion, the function of detecting complex structural variation types is achieved. In the genotyping part, a genotyping method based on distance similarity clustering and sequence similarity classification is proposed. By combining distance clustering with hierarchical clustering and using the minimum edit distance, haplotypes are divided more finely to solve the accuracy and efficiency problems of structural variation typing. The present invention provides a complete set of processes and methods for end-to-end splicing genome detection of structural variation, which can accurately and efficiently detect common structural variation and complex structural variation types of the genome and accurately and efficiently complete the genotyping function, filling the gaps and deficiencies of existing methods in end-to-end genome detection of structural variation; at the same time, the clustering and filtering methods proposed in the present invention can effectively reduce false positives, avoid the influence of noise signals, and improve detection accuracy; the complex variation structural variation detection method proposed can quickly realize the identification of complex structural variation, and provide clues to reveal the manifestations and laws of complex variation patterns; the proposed genotyping method effectively solves the accuracy and efficiency problems of structural variation typing. The overall solution is fast, efficient, accurate, and has great application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 It is the overall flow chart of the present invention;

[0017] Figure 2 This is the flow chart of mutation signal clustering;

[0018] Figure 3 Flow chart of genotyping. DETAILED DESCRIPTION

[0019] Specific implementation method 1: The specific process of the variation detection method based on end-to-end genome assembly in this implementation method is as follows:

[0020] In order to obtain the final variation file from the spliced ​​genome alignment file, the algorithm mainly includes three steps: variation signal processing part, complex structural variation detection part and genotyping part.

[0021] In order to obtain the bam file for analysis, it is necessary to first obtain the quality-controlled sequencing data in fastq format from high-throughput sequencing platforms such as Illumina, PacBio, etc., and then assemble and type them through splicing tools such as Hifiasm, Peregrine, Wtdbg2, etc. to form two haploid genome splicing fasta files, and then use alignment tools such as minimap2, NGMLR, Winnowmap and LRA to align the assembled genome with the reference genome to generate a SAM format file, and then use samtools to convert the SAM into a binary BAM file.

[0022] In the variation signal processing, the overall process is to obtain all the comparison information records from the bam file, divide the information in the comparison record into intra signals and inter signals according to their flag bits, and perform signal processing on them respectively to obtain five basic variation signals.

[0023] To realize the variation signal processing, it is divided into three steps: comparison information filtering, variation signal extraction and variation signal clustering.

[0024] Step 1: Process the variant signal to obtain the processed variant signal; the specific process is:

[0025] Step 11, obtain BAM format comparison information;

[0026] Step 12: Filter the BAM format comparison information obtained in step 11 to obtain filtered BAM format comparison information;

[0027] Step 13: extract the variation signal from the filtered BAM format comparison information;

[0028] Step 14: Cluster the extracted variant signals and output the clustering results;

[0029] Step 15: Filter all the output cluster results to obtain the final cluster result, that is, integrate the signals contained in each cluster into a structural variation signal for output;

[0030] Step 2: perform mixed sorting on the processed variation signals obtained in step 1, perform spatial distribution clustering on the sorted variation signals, and determine whether the variation signals after spatial distribution clustering are potential complex variations. If so, it is determined to be a complex variation signal, and a complex variation signal is obtained, and step 3 is executed; if not, it is determined to be a simple variation signal, and step 3 is executed. A simple variation signal is a variation signal of insertion, deletion, inversion, translocation, and duplication;

[0031] Step 3: Perform genotyping on the variant signal obtained in step 2.

[0032] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the BAM format comparison information is obtained in step 1; the specific process is:

[0033] 1) Obtain sequencing data files in fastq format from a high-throughput sequencing platform;

[0034] Use Illumina and PacBio to obtain fastq format sequencing data files from high-throughput sequencing platforms;

[0035] 2) Use a splicing tool to process the sequencing data file in fastq format to obtain two haplotype assembled genomes after typing. The format of the two haplotype assembled genomes is fasta format;

[0036] Typing refers to the genetic composition of diploid organisms, usually one gene comes from the father and one from the mother, and it is necessary to determine whether the gene comes from the father or the mother;

[0037] Use assembly and typing tools such as Hifiasm, Peregrine, Wtdbg2, etc. to process the sequencing data files in fastq format and obtain two haploid genome splicing fasta files;

[0038] 3) Use the alignment tool to align the two haplotype assembled genomes with the reference genome, and generate two SAM format alignment information files;

[0039] Align the assembled genome with the reference genome using alignment tools such as minimap2, NGMLR, Winnowmap, and LRA to generate SAM format files;

[0040] 4) Convert the two SAM format comparison information files into binary BAM format comparison information files respectively.

[0041] Use samtools to convert SAM format comparison information into binary BAM format comparison information.

[0042] The other steps and parameters are the same as those in the first embodiment.

[0043] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that in step one or two, the BAM format comparison information obtained in step one or two is filtered to obtain the filtered BAM format comparison information; the specific process is:

[0044] Classify the comparison information record according to the flag bit in the BAM format comparison information to determine the type of the comparison information record;

[0045] The bam format comparison information consists of many lines, each line represents a record;

[0046] The types of comparison information records can be divided into primary comparison, secondary comparison or supplementary comparison;

[0047] Each alignment information records the read / contig name from which it comes, the name of the aligned chromosome, the starting position aligned to the reference genome, the alignment quality, the aligned CIGAR string, the aligned sequence, the alignment supplementary information, etc.

[0048] The supplementary information contains other tags to provide additional information. For example, the SA tag records the alignment information of all the splits that make up the same spliced ​​sequence fragments. That is, through the SA tag of the main alignment, according to the same contig name (the first column), the alignment information record of all the fragments that make up the same contig is obtained, as given in Formula 1:

[0049] contig={seg 1 ,seg 2 ,…,seg j ,…,seg n}, 1 <i≤n (1)

[0050] The fasta file is input into the alignment software to obtain the bam file. The alignment software may split the fasta records with the same name into several alignment records. These several alignment records come from the same fasta record, that is, from the same contig assembly fragment. Alignment means finding the position of this gene sequence in the reference genome. To clearly represent this position, it is necessary to record which chromosome (chromosome name) this gene sequence is aligned to and which position (chromosome start position) on the chromosome. To indicate the credibility of this alignment, there will be an alignment quality to measure. The higher the quality, the better and more credible the alignment. The CIGAR string of the alignment indicates the form of this alignment, that is, through what operations this aligned sequence will be restored to the reference genome sequence, such as operations like adding (insertion, I), deleting (deletion, D), or being identical (M), etc. The aligned sequence refers to the specific sequence of the gene composition of this alignment, such as AGCT. The gene sequence is formed by the permutation and combination of the four characters AGCT.

[0051] Before performing conventional variant signal processing, usually there is a simple filtering of the alignment quality of the alignment fragment information, screening out the alignment fragment information with an alignment quality greater than a certain threshold (usually 20). However, due to the phenomenon of redundant alignment in the alignment information, that is, the same gene sequence may be aligned to multiple positions in the genome to generate multiple alignment information records. Therefore, here a filtering condition between alignment information records is provided to obtain maximally non-overlapping alignment fragments to reduce false positives.

[0052] The filtering conditions are:

[0053]

[0054] Among them, contig is the alignment information of the assembly fragment, and seg i is the alignment information record of the i-th fragment sorted according to the relative start position of the fragment in the contig. ReadEnd segi represents the end position of seg i in the contig. represents the end position of seg i-1 in the contig; 1 < i ≤ n, where n represents the total number of fragments in the contig.

[0055] Select the alignment information as the filtered BAM format alignment information.

[0056] Other steps and parameters are the same as those in the specific implementation method one or two.

[0057] Specific implementation method 4: This implementation method is different from one of the specific implementation methods 1 to 3 in that in step 13, the variation signal in the filtered BAM format comparison information is extracted; the specific process is:

[0058] The comparison information is divided into intra signal and inter signal according to the flag bit in the filtered BAM format comparison information;

[0059] Perform signal processing on the intra signal to obtain two types of variation signals: insertion and deletion;

[0060] Perform signal processing on the inter signal to obtain five types of variation signals: insertion, deletion, inversion, translocation, and duplication;

[0061] The intra signal is processed to obtain two kinds of variation signals, insertion and deletion. The specific process is as follows:

[0062] If the filtered seg i If the number before the I character in the CIGAR character position in is greater than 30bp, it is considered that there is an insertion signal, and the chromosome name where the insertion signal is located, the starting position of the insertion signal on the reference genome, the size of the insertion signal, and the insertion signal sequence are recorded as an insertion signal; otherwise, there is no insertion signal;

[0063] If the filtered seg i If the number before the D character in the CIGAR character position in is greater than 30bp, it is considered that there is a deletion signal, and the chromosome name where the deletion signal is located, the starting position of the deletion signal on the reference genome, the ending position of the deletion signal on the reference genome, and the size of the deletion signal are recorded as a deletion signal; otherwise, there is no deletion signal.

[0064] In the specific signal extraction part, the intra signal refers to the signal directly obtained from each seg after filtering above. i Insertion and deletion signals greater than 30bp are obtained from the CIGAR character position in the alignment information record. The format of CIGAR is a number plus type, such as "87H23M2I38M60I29M39D200H". If the number before the I character in the CIGAR character position is greater than 30bp, it is considered that there is an insertion signal, and its chromosome name, starting position on the reference genome, insertion size, insertion sequence, etc. are recorded as an insertion signal. For example, "60I" indicates that there is an insertion signal with a length of 60 here, and the rest of the information can be directly obtained from the alignment record; if the number before the D character in the CIGAR character position is greater than 30bp, it is considered that there is a deletion signal, and its chromosome name, starting position on the reference genome, ending position on the reference genome, deletion size, etc. are recorded as a deletion signal.

[0065] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0066] Specific implementation method 5: This implementation method is different from specific implementation methods 1 to 4 in that the inter signal is processed to obtain five types of variation signals: insertion, deletion, inversion, translocation, and duplication; the specific process is:

[0067] Inter signals are implemented by setting heuristic rules based on the difference in the front-to-back distances between the filtered segs on the spliced ​​fragment contig, as well as the difference in chromosome numbers and physical distances compared on the reference genome. For example, for the five basic variation signals of different types (insertion, deletion, inversion, translocation, and duplication), different distance threshold judgment criteria are formulated based on experience and previous data analysis results: if the front-to-back distance of a signal on the spliced ​​fragment is within a specific range, and it is compared to a specific chromosome on the reference genome and the distance meets another condition, it is considered to meet the extraction requirements and the extraction process is triggered. Specifically, it is through the read dis and ref dis And the corresponding information extracts five specific mutation signals, as follows:

[0068] Traverse the two adjacent segments that make up the contig, denoted as seg 1 and seg 2 , where the fragment with the smaller starting coordinate is designated as the first aligned fragment;

[0069]

[0070] in, Indicates seg 1 The starting position on the spliced ​​fragment, Indicates seg 2 The starting position on the spliced ​​fragment;

[0071] Indicates seg 1 The end position on the spliced ​​fragment, Indicates seg 2 The termination position on the spliced ​​fragment;

[0072] Indicates seg 1 The starting position on the reference genome, Indicates seg 2 The starting position on the reference genome;

[0073] Indicates seg 1The stop position on the reference genome, Indicates seg 2 The stop position on the reference genome;

[0074] read dis Indicates seg 2 The starting position on the spliced ​​fragment is the same as seg 1 The distance between the end positions on the spliced ​​fragments, ref dis Indicates seg 2 The starting position on the reference genome is the same as seg 1 The distance between the stop positions on the reference genome, OLP ref Indicates seg 1 The end position on the reference genome is the same as seg 2 The overlap size of the starting position on the reference genome on the reference genome;

[0075] 1) For the insertion (INS) mutation signal, if the insertion (INS) mutation signal meets the following preconditions:

[0076]

[0077] in, Indicates seg 1 Chromosome name, Indicates seg 2 The chromosome name in formula (4) indicates that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are consistent;

[0078] And satisfy:

[0079]

[0080] Then extract the insert (INS) signal:

[0081]

[0082] Among them, Chrom represents the chromosome name identifier where the extracted signal is located. INS Indicates the chromosome name where the extracted insertion signal is located (the chromosome name recorded in the first alignment information or the second alignment information), Start indicates the starting position of the extracted signal, RefEnd seg1 Indicates seg 1 The end position of the reference genome, RefEnd seg2 Indicates seg 2The end position of the reference genome in , SVlen represents the variant length, Seg represents the sequence information identifier of the insertion signal, Seq INS Indicates the sequence corresponding to the insertion signal;

[0083] If formulas (4) and (5) are not satisfied, the insertion (INS) signal is not extracted;

[0084] 2) For the deletion (DEL) mutation signal, if the deletion (DEL) mutation signal meets the following preconditions:

[0085]

[0086] in, Indicates seg 1 Chromosome name, Indicates seg 2 The chromosome name in formula (4) indicates that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are consistent;

[0087] And satisfy:

[0088]

[0089] Then extract the delete (DEL) signal:

[0090]

[0091] Among them, Chrome DEL Indicates the chromosome name where the extracted deletion signal is located (the chromosome name recorded in the first alignment information or the second alignment information), End indicates the end position mark of the extracted signal, and other marks are defined as above;

[0092] If formulas (4) and (7) are not satisfied, the delete (DEL) signal is not extracted;

[0093] 3) For the duplication (DUP) variant signal, if the duplication (DUP) variant signal meets the following preconditions:

[0094]

[0095] in, Indicates seg 1 Chromosome name, Indicates seg 2 The chromosome name in formula (4) indicates that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are consistent;

[0096] And satisfy:

[0097]

[0098] Then extract the duplicate (DUP) signal:

[0099]

[0100] Among them, ReadEnd seg1 Indicates seg 1 The end position of the spliced ​​fragment, ReadStart seg2 Indicates seg 2 The starting position of the spliced ​​fragment, RefStrart seg2 Indicates seg 2 The starting position of the reference genome in Chrom DUP Indicates the chromosome name where the extracted DUP signal is located (the chromosome name recorded in the first alignment information or the second alignment information), Th olp Indicates the distance threshold for extracting repeated DUP signals (the default value is 5000bp). Different distance thresholds are set according to different types of mutations to more accurately identify each mutation;

[0101] If formulas (4) and (9) are not satisfied, the repeated (DUP) signal is not extracted;

[0102] 4) For inversion (INV) mutation signals;

[0103] 41) If the inversion (INV) mutation signal meets the following preconditions:

[0104]

[0105] in, Indicates seg 1 Chromosome name, Indicates seg 2 The chromosome name in formula (4) indicates that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are consistent;

[0106] And satisfy:

[0107] seg 1 and seg 2 The comparison direction in the record is different, seg 1 and seg 2 There is a reverse alignment and a forward alignment. The seg of the forward alignment is recorded as fwd, and the seg of the reverse alignment is recorded as rev;

[0108] |ReadEnd rev -ReadStart fwd |≤100bp (11)

[0109] Then extract the inversion (INV) signal:

[0110]

[0111] If formulas (4) and (11) are not satisfied, the inversion (INV) signal is not extracted;

[0112] 42) If the inversion (INV) mutation signal meets the following preconditions:

[0113]

[0114] in, Indicates seg 1 Chromosome name, Indicates seg 2 The chromosome name in formula (4) indicates that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are consistent;

[0115] And satisfy:

[0116] seg 1 and seg 2 The comparison direction in the record is different, seg 1 and seg 2 There is a reverse alignment and a forward alignment. The seg of the forward alignment is recorded as fwd, and the seg of the reverse alignment is recorded as rev;

[0117] |ReadEnd fwd -ReadStart rev |≤100bp (13)

[0118] Extract the INV signal:

[0119]

[0120] If formulas (4) and (13) are not satisfied, the inversion (INV) signal is not extracted;

[0121] Among them, ReadEnd rev ReadStart is the end position of the spliced ​​fragment in the reverse alignment information. fwd RefStart is the starting position of the spliced ​​fragment in the forward alignment information. rev RefStart is the starting position of the reference genome in the reverse alignment information. fwdis the starting position of the reference genome in the forward alignment information, ReadEnd fwd ReadStart is the end position of the spliced ​​fragment in the forward alignment information. rev RefEnd is the starting position of the spliced ​​fragment in the reverse alignment information, fw d is the reference genome end position in the forward alignment information, RefEnd rev The reference genome end position in the reverse alignment information; Chrom INV Start indicates the chromosome name where the INV mutation signal occurs; INV Indicates the start position of the INV mutation signal, End INV The mark indicating the end position of the INV variation signal;

[0122] 5) For translocation (TRA) mutation signal;

[0123] 51) If the translocation (TRA) mutation signal meets the following preconditions:

[0124]

[0125] Among them, formula (15) shows that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are inconsistent;

[0126] And satisfy:

[0127] |ReadEnd seg1 -ReadStart seg2 |≤1000bp (16)

[0128] Then extract the TRA signal:

[0129]

[0130] If formulas (15) and (16) are not satisfied, the translocation (TRA) signal is not extracted;

[0131] 52) If the translocation (TRA) mutation signal meets the following preconditions:

[0132]

[0133] Among them, formula (15) shows that seg 1 The chromosome names recorded in seg 2 The chromosome names recorded are inconsistent;

[0134] And satisfy:

[0135] |ReadEnd seg2-ReadStart seg1 |≤1000bp (18)

[0136] Then extract the translocation (TRA) signal:

[0137]

[0138] If formulas (15) and (18) are not satisfied, the translocation (TRA) signal is not extracted;

[0139] Among them, Chrome 1 The identifier of the chromosome where the TRA mutation signal occurs, Pos 1 Indicates the starting position of the TRA mutation signal. 2 The identifier of the termination chromosome name where the TRA mutation signal occurs, Pos 2 The marker indicates the end position where the TRA mutation signal occurs.

[0140] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0141] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the extracted variant signals are clustered in step 14 to output the clustering results; the specific process is:

[0142] The existence of mutation signal clustering is necessary. There are two reasons: First, due to the deviation of the genome assembly and genome sequence alignment software and the existence of a large number of repeated regions in the genome itself, there is a section containing multiple mutations of the same type, and most of them may be false positives; second, there are some mutations of the same type that are distributed close to each other and can be aggregated into one mutation. Here, an unsupervised clustering and filtering method based on feature similarity is provided: for each mutation type, mutations with close sites and similar features are clustered together based on the distance distribution to form an average mutation to achieve the filtering effect. For abnormally prominent clustered sites of mutations of the same type, they are discarded because the credibility is not high under the current technical flow, which increases the probability of false positives.

[0143] (1) For the extracted deletion (DEL) signal, insertion (INS) signal, duplication (DUP) signal, inversion (INV) signal, and translocation (TRA) signal;

[0144] Initialize the deletion (DEL) cluster container; initialize the insertion (INS) cluster container; initialize the duplication (DUP) cluster container; initialize the inversion (INV) cluster container; initialize the translocation (TRA) cluster container;

[0145] Initialize clustering containers, and create corresponding clustering containers according to the input mutation feature set and the predetermined mutation types (such as DEL, INS, DUP, INV, TRA);

[0146] (2) Sort the variant signals in each variant type set in ascending order according to the starting position where they appear in the reference genome;

[0147] Sorting the signal set of each variant type according to a preset sorting rule so that the variants are sorted in ascending order according to the starting position where they appear in the reference genome;

[0148] (3) For the signals in the mutation signal set of the deletion (DEL) mutation type, perform the following operations in sequence:

[0149] a. Initialize an empty cluster and add the first signal in the mutation signal set of the deletion mutation type as the starting signal, and execute d;

[0150] b. Calculate the similarity between the current signal and the last signal in the current cluster;

[0151] c. If the similarity is less than a predetermined threshold, the current cluster is added to the delete (DEL) cluster container and the current signal is added to the new cluster;

[0152] If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster;

[0153] d. Move to the next signal and repeat bc until all signals in the mutation signal set of the deletion (DEL) mutation type are determined;

[0154] (4) For the signals in the mutation signal set of the insertion (INS) mutation type, perform the following operations in sequence:

[0155] a. Initialize an empty cluster and add the first signal in the mutation signal set of the inserted mutation type as the starting signal, and execute d;

[0156] b. Calculate the similarity between the current signal and the last signal of the current cluster;

[0157] c. If the similarity is less than a predetermined threshold, the current cluster is added to the insert (INS) cluster container and the current signal is added to the new cluster;

[0158] If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster;

[0159] d. Move to the next signal and repeat bc until all signals in the mutation signal set of the insertion (INS) mutation type are determined;

[0160] (5) For the signals in the mutation signal set of the duplication (DUP) mutation type, perform the following operations in sequence:

[0161] a. Initialize an empty cluster and add the first signal in the set of mutation signals of the repeated mutation type as the starting signal, and execute d;

[0162] b. Calculate the similarity between the current signal and the last signal of the current cluster;

[0163] c. If the similarity is less than a predetermined threshold, the current cluster is added to the duplicate (DUP) cluster container and the current signal is added to the new cluster;

[0164] If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster;

[0165] d. Move to the next signal and repeat bc until all signals in the set of variant signals of the duplication (DUP) variant type are determined;

[0166] (6) For the signals in the mutation signal set of the inversion (INV) mutation type, the following operations are performed in sequence:

[0167] a. Initialize an empty cluster, add the first signal in the mutation signal set of the inversion mutation type as the starting signal, and execute d;

[0168] b. Calculate the similarity between the current signal and the last signal of the current cluster;

[0169] c. If the similarity is less than a predetermined threshold, the current cluster is added to the inversion (INV) cluster container, and the current signal is added to the new cluster;

[0170] If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster;

[0171] d. Move to the next signal and repeat bc until all signals in the variant signal set of the inversion (INV) variant type are determined;

[0172] (7) For the signals in the mutation signal set of the translocation (TRA) mutation type, the following operations are performed in sequence:

[0173] a. Initialize an empty cluster, add the first signal in the set of mutation signals of the translocation mutation type as the starting signal, and execute d;

[0174] b. Calculate the similarity between the current signal and the last signal of the current cluster;

[0175] c. If the similarity is less than a predetermined threshold, the current cluster is added to the translocation (TRA) cluster container, and the current signal is added to the new cluster;

[0176] If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster;

[0177] d. Move to the next signal and repeat bc until all signals in the set of variant signals of the translocation (TRA) variant type are determined;

[0178] (8) Obtain all clustering results of all variant type signals.

[0179] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0180] Specific implementation method 7: This implementation method is different from any one of specific implementation methods 1 to 6 in that the similarity between each signal and the last signal of the current cluster is calculated in b; the specific process is:

[0181] 1. For insertion, deletion, inversion, and duplication mutation signals:

[0182] 1) Calculate the location similarity S pos Scoring; the specific process is:

[0183] Position similarity S pos The score calculation formula is as follows:

[0184] S pos =|Start 1 -Start 2 | (20)

[0185] Among them, Start 1 、Start 2 They are respectively the starting positions of the signals in the two variant signals. The position difference is extremely critical in the similarity assessment.

[0186] 2) Calculate the similarity score S of the mutation signal size size ; The specific process is:

[0187] The mutation lengths of the two mutation signals are SVlen 1 and SVlen 2 , the variation length size similarity score S of the variation signal size Calculated by the following formula:

[0188]

[0189] For inversion and duplication mutation signals, the SVlen can be calculated by the difference between the mutation end position and the mutation start position.

[0190] This simple proportional formula can intuitively reflect the degree of fit between the sizes of two variants. The value range is between 0 and 1. The closer it is to 1, the closer the sizes of the two variants are.

[0191] 3) Calculate the chromosome similarity score S Chr ; The specific process is:

[0192] Assume that the chromosome numbers of the two mutation signals are Chrom 1 and Chrome 2 , chromosome similarity score S Chr Calculated by the following discrete function:

[0193]

[0194] When the chromosome names of the two variant signals are the same, the chromosome similarity score is 1;

[0195] When the chromosome names of the two variant signals are different, the chromosome similarity score is 0;

[0196] 4) The position similarity S pos Score, variation size similarity score S size Score, chromosome similarity score S Chr Perform weighted summation to obtain a comprehensive similarity score S; expressed as:

[0197] S=S Chr ×(w size ×S size +w pos ×S pos ) (twenty three)

[0198] Among them, w pos 、w size are the weight coefficients of the position similarity score and the variation size similarity score, respectively, and w pos +w size =1, for example, 0.4, 0.6. The weight can be adjusted according to the specific analysis scenario and data characteristics to make the formula better adapt to actual needs and accurately measure the similarity between variant features.

[0199] 2. For translocation mutation signals:

[0200]

[0201] Among them, S Chr1 represents the starting chromosome similarity score, Indicates the name of the starting chromosome of the first signal, indicates the name of the starting chromosome of the second signal;

[0202]

[0203] Among them, S Chr2 represents the termination chromosome similarity score, The name of the chromosome that terminates the first signal, The name of the chromosome that terminates the second signal;

[0204] The starting chromosome similarity score S Chr1 , end chromosome similarity score S Chr2 Perform weighted summation to obtain a comprehensive similarity score S; expressed as:

[0205]

[0206] in, refers to the starting position of the second signal, Refers to the starting position of the first signal, Refers to the end position of the second signal Refers to the end position of the first signal;

[0207] w pos1 is the weight coefficient of the starting position, w pos2 are the weight coefficients of the end position, and the default values ​​are 0.5 and 0.5 respectively.

[0208] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0209] Specific implementation eight: This implementation differs from specific implementations one to seven in that in step one-five, all the output cluster results are filtered to obtain the final cluster result, that is, the signals contained in each cluster are integrated into a structural variation signal for output; the specific process is:

[0210] 1) Set the maximum number of clusters max_cluster, the default value is 5;

[0211] 2) If the number of clusters in a cluster is greater than max_cluster, the corresponding cluster is discarded; 3)

[0213] 31) Clustering results for insertion mutation types:

[0214] For each clustering cluster of the clustering result of the inserted mutation type, if the number of clusters in the clustering cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and average length of all mutation data in the clustering cluster; the specific process is:

[0215] The mutation start positions of each mutation signal in the δth cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ = 1, 2, ..., N, N represents the total number of clusters in the clustering results of the insertion mutation type;

[0216] The lengths of all variants in the δth cluster are summed up and divided by the number of variants to calculate the average length;

[0217] Keep the other information in the δth cluster the same as the first mutation signal in the cluster;

[0218] The other information is the chromosome name, the name from the contig, and the inserted sequence;

[0219] The average mutation start position, average length and other information are output as the δth clustering result;

[0220] Until the output results of all clustering clusters of the clustering results of the insertion mutation type are obtained, that is, the final insertion mutation signal;

[0221] 32) Clustering results for deletion mutation types:

[0222] For each clustering cluster of the clustering results of the deleted variant type, if the number of clusters in the clustering cluster is less than or equal to max_cluster, calculate the average variant signal starting position, average variant signal ending position and average length of all variant data in the clustering cluster; the specific process is:

[0223] The mutation start positions of each mutation signal in the δ′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ′=1,2,…,N′, N′ represents the total number of clusters in the clustering result of deleting the mutation type;

[0224] The mutation termination positions of each mutation signal in the δ′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position;

[0225] The lengths of all variants in the δ′th cluster are summed up and divided by the number of variants to calculate the average length.

[0226] Keep the other information in the δ′th cluster the same as the first mutation signal in the cluster;

[0227] The other information is the chromosome name and the name from the contig;

[0228] The average mutation start position, average mutation end position, average length and other information are output as the δ′th cluster;

[0229] Until the output results of all clustering clusters of the clustering results of the deleted mutation type are obtained, that is, the final deleted mutation signal;

[0230] 33) Clustering results for inversion mutation types:

[0231] For each cluster of the inversion mutation type clustering result, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and the average mutation signal ending position of all mutation data in the cluster; the specific process is:

[0232] The mutation start positions of each mutation signal in the δ″th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ″=1,2,…,N″, N″ represents the total number of clusters in the clustering results of the inversion mutation type;

[0233] The mutation termination positions of each mutation signal in the δ″th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position;

[0234] Keep the other information in the δ″th cluster the same as the first mutation signal in the cluster;

[0235] The other information is the chromosome name and the name from the contig;

[0236] The average mutation start position, average mutation end position and other information are used as the output result of the δ″th cluster;

[0237] Until the output results of all clustering clusters of the clustering results of the inversion mutation type are obtained, that is, the final inversion mutation signal;

[0238] 34) Clustering results for translocation variant types:

[0239] For each cluster of the clustering results of the translocation variant type, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and the average mutation signal ending position of all mutation data in the cluster; the specific process is:

[0240] The mutation start positions of each mutation signal in the δ″′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ″′=1,2,…,N″′, N″′ represents the total number of clusters in the clustering results of the translocation mutation type;

[0241] The mutation termination positions of each mutation signal in the δ″′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position;

[0242] Keep the other information in the δ″′th cluster the same as the first mutation signal in the cluster;

[0243] The other information is the name of the starting chromosome, the name of the ending chromosome, and the name from the contig;

[0244] The average mutation start position, average mutation end position and other information are output as the δ″′th clustering result;

[0245] Until the output results of all clustering clusters of the clustering results of the translocation variation type are obtained, that is, the final translocation variation signal;

[0246] 35) Clustering results for repeated variation types:

[0247] For each cluster of the clustering results of the repeated mutation type, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and the average mutation signal ending position of all mutation data in the cluster; the specific process is:

[0248] The first The mutation start positions of each mutation signal in each cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; N″′ represents the total number of clusters in the clustering results of repeated variation types;

[0249] The first The mutation termination positions of each mutation signal in each cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position;

[0250] The first The other information in the cluster remains the same as the first mutation signal in the cluster;

[0251] The other information is the chromosome name and the name from the contig;

[0252] The average mutation start position, average mutation end position and other information are used as the Cluster output results;

[0253] Until the output results of all clustering clusters of the clustering results of the repeated variation type are obtained, that is, the final repeated variation signal.

[0254] Because some of the signals are numerical types, such as the chromosome starting position and length, which can be averaged in the cluster, but some information, such as the chromosome name, is string type and cannot be averaged, so you only need to keep the same content as the first mutation signal in the cluster, because in the same cluster, these non-numerical types are consistent;

[0255] Filtering mainly includes the following steps:

[0256] 1) In the subsequent processing stage of cluster analysis, strict quantity screening operations must be performed for each generated cluster. The default maximum number of clusters is 5, which is used as the key judgment criterion. After the clustering operation is completed and a series of clusters are obtained, the number of clusters is counted and judged.

[0257] 2) If the number of clusters obtained by statistics exceeds the maximum number of clusters, such clusters will be directly discarded and will no longer participate in the subsequent result generation process. The purpose of this is to avoid too many and too complex clusters introducing unnecessary noise data, interfering with the final analysis results, and ensuring the simplicity and reliability of the output results.

[0258] 3) On the contrary, if the number of clusters is less than the default number of clusters, a detailed result extraction process is started. At this time, the average breakpoint position and average length are calculated for all variant data in the cluster. Specifically, the breakpoint positions (breakpoints are the starting position of chromosomes or the ending position of chromosomes) of each variant in the cluster are added up one by one, and then divided by the number of variants to accurately obtain the average breakpoint position; similarly, the lengths of all variants are added up and divided by the number of variants to calculate the average length. Finally, the two calculated average breakpoint positions and average lengths and other information that remains the same as the first signal are used as the final output results corresponding to the cluster.

[0259] 4) Repeat 2) to 3) until all clusters are filtered.

[0260] In clustering, similarity calculation is based on a predefined feature similarity formula, which calculates similarity scores based on the attributes of the features (such as position, size, chromosome); the threshold is a user-adjustable parameter used to define the similarity judgment criteria between features. In actual operation, when the similarity score between two features calculated based on the similarity formula is higher than the threshold, the two features are judged to be similar and classified into the same cluster; conversely, if the score is lower than the threshold, the two features are judged to be very different and are divided into different clusters. This method can automatically adapt to feature sets of different sizes and has low computational complexity.

[0261] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0262] Specific implementation method nine: This implementation method is different from specific implementation methods one to eight in that, in step two, the processed variation signals obtained in step one are mixed and sorted, and the sorted variation signals are spatially distributed and clustered to determine whether the variation signals after spatial distribution and clustering are potential complex variations. If so, it is determined to be a complex variation signal, and a complex variation signal is obtained, and step three is executed; if not, it is determined to be a simple variation signal, and step three is executed. A simple variation signal is a variation signal of insertion, deletion, inversion, translocation, and duplication. The specific process is:

[0263] For the complex variation part, the basic variation signals obtained in the previous stage are mixed and sorted according to the position of the variation signals on the reference genome to ensure that the subsequent clustering analysis can effectively group based on spatial position; then clustering is performed according to distance, and complex structural variation fusion is further performed for potential complex variation sites that meet the corresponding pattern to obtain the final complex variation form.

[0264] 1) Use merge sort or quick sort to sort the processed variant signals obtained in step 1 in ascending order according to their starting positions on the reference genome;

[0265] After the first step, there is a container for each variant type. Now we need to mix these types together and re-sort them according to their positions in the reference genome. The method can be merge or quick sort, but the result after sorting is certain, that is, all variants are sorted in ascending order according to their starting position on the reference genome;

[0266] The interval distance parameter Dist is introduced and set to 1000bp by default. This parameter serves as a key scale indicator for subsequent determination of potential complex variation regions.

[0267] 2) Taking the genome position where the first mutation signal after sorting is located as the center, extend Dist / 2 base pairs forward from the center, and extend Dist / 2 base pairs backward from the center to form a Dist scanning window of length;

[0268] Acquire all mutation signals within the scanning window;

[0269] Starting from the first variant signal after sorting, with the genomic position of the variant signal as the center, a genomic interval of length Dist is constructed as the scanning window. Specifically, the window extends Dist / 2 base pairs forward and the same number of base pairs backward. Subsequently, all variant signals falling within the scanning window are checked in order.

[0270] 3) Determine the mutation type to which all mutation signals in the scanning window belong;

[0271] The types of variants include insertion, deletion, inversion, translocation, or duplication;

[0272] When the number of variant signal types in a certain scanning window exceeds two, the scanning window is identified as a potential complex variant region; otherwise, the scanning window is not a potential complex variant region;

[0273] For the variation signals covered by each scanning window, the number of different types of structural variation signals is strictly counted, including insertion, deletion, inversion, translocation and duplication. When the number of different types of structural variation signals counted in a scanning window exceeds two, this region is identified as a potential complex variation region.

[0274] 4) Extract detailed information of each complex variation signal in the potential complex variation region. The detailed information includes the variation type, variation chromosome name, and variation start position of each variation signal. The format of the complex variation signal is (CSV, chromosome name, variation type, variation start position);

[0275] "CSV" is a universal comma-separated value format identifier, which facilitates data storage, transmission and subsequent analysis and processing.

[0276] Once a certain interval is determined to be a potential complex variation region, the complex variation pattern fusion process is immediately started. The system comprehensively collects detailed information on all variation signals in the region, covering key elements such as variation type, chromosome name, and starting position, and integrates and generates a set of structured complex variation description information. The format is defined as (CSV, variation type, variation starting position), where "CSV" is a universal comma-separated value format identifier, which is convenient for data storage, transmission, and subsequent analysis and processing.

[0277] 5) Move the scanning window back one variant signal in order, and repeat 1)-4) until all complex variant signals after mixed sorting are traversed and processed.

[0278] After completing the scanning, judgment, fusion and output process of the window corresponding to the first variant signal, the scanning window is moved back one variant signal in sequence, and the above whole process is repeated seamlessly until all the sorted variant signals are traversed and processed.

[0279] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0280] Specific implementation method ten: This implementation method is different from the specific implementation methods one to nine in that in step three, the variant signal obtained in step two is subjected to genotyping; the specific process is as follows:

[0281] Genotyping is a key step in genomic research, which is used to analyze the distribution of genotypes in a population and their relationship with phenotypes. However, due to the complexity of structural variation, traditional genotyping methods are difficult to accurately handle. Existing technologies rely on manual annotation or simple rules, resulting in low typing accuracy, especially when processing high-throughput data, the efficiency is significantly reduced. The present invention provides a genotyping method based on distance similarity clustering and sequence similarity classification, which solves the accuracy and efficiency problems of complex variation typing by combining distance clustering with hierarchical clustering. The method includes two main steps: first, distance clustering of variant sites and variant lengths according to the variant type; second, hierarchical clustering based on sequence similarity is performed within each distance clustering group, and the haplotype determination result is calculated in combination with the minimum edit distance. The specific steps are as follows:

[0282] Step 31:

[0283] All the variation signal information obtained after the above steps includes the chromosome name, variation start position, variation end position, variation length and sequence information, etc. The specific information varies depending on the type of variation signal.

[0284] Distance clustering is performed on different types of variation signals obtained in step 2; the specific process is as follows:

[0285] For each mutation type, add the initial mutation signal to the clustering cluster, and traverse the next mutation signal sequentially until the corresponding type of mutation signal completes the clustering:

[0286] For insertion (INS) type mutation signals, three conditions are checked: first, the chromosome name where the current mutation signal is located is the same as the chromosome name where the mutation is located in the current cluster; second, the Euclidean distance between the chromosome starting position in the current mutation signal and the chromosome starting position of the last mutation signal in the current cluster is less than 200bp; third, the difference between the mutation length of the current mutation signal and the mutation length of the last mutation signal in the current cluster is less than 50bp; if all three conditions are met, the current mutation signal is added to the current cluster; otherwise, a new cluster is opened and the current mutation signal is added to the new cluster;

[0287] For mutation signals of deletion (DEL), inversion (INV) or duplication (DUP) types, three conditions are checked: first, the chromosome name where the current mutation signal is located is the same as the chromosome name where the mutation is located in the current cluster; second, the Euclidean distance between the mutation start position of the current mutation signal and the mutation start position of the last mutation signal in the current cluster is less than 200bp; third, the Euclidean distance between the mutation end position in the current mutation signal and the mutation end position in the last mutation signal in the current cluster is less than 200bp; if all three conditions are met, the current mutation signal is added to the current cluster; otherwise, a new cluster is opened and the current mutation signal is added to the new cluster;

[0288] For the translocation (TRA) type mutation signal, four conditions are checked: first, the chromosome name where the current mutation signal starts is the same as the chromosome name where the last mutation signal starts in the current cluster; second, the chromosome name where the current mutation signal ends is the same as the chromosome name where the last mutation signal ends in the current cluster; third, the Euclidean distance between the mutation start position of the current mutation signal and the mutation start position of the last mutation signal in the current cluster is less than 200bp; fourth, the Euclidean distance between the mutation end position of the current mutation signal and the mutation end position of the last mutation signal in the current cluster is less than 200bp; if all four conditions are met, the current mutation signal is added to the current cluster; otherwise, a new cluster is opened and the current mutation signal is added to the new cluster;

[0289] For the mutation signal of complex mutation (CSV), three conditions are checked: first, the chromosome name where the current mutation signal starts is the same as the chromosome name where the last mutation starts in the current cluster; second, the Euclidean distance between the mutation start position of the current mutation signal and the last mutation start position in the current cluster is less than 1000bp; third, the mutation type of the current mutation signal is the same as the mutation type of the last mutation signal in the current cluster; if all three conditions are met, the current signal is added to the current cluster; otherwise, a new cluster is opened and the current signal is added to the new cluster;

[0290] Through the above steps, each cluster is used as a distance cluster indicating potential identical haplotype variations;

[0291] Step 32

[0292] If the signal within the distance clustering cluster is an insertion mutation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the insertion mutation sequences of any two mutation signals;

[0293] If the signal within the distance clustering cluster is a deletion mutation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the deletion mutation sequences of any two mutation signals;

[0294] If the signal within the distance clustering cluster is a repeated variation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the repeated variation sequences of any two variation signals;

[0295] If the signal within the distance clustering cluster is an inversion mutation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the inversion mutation sequences of any two mutation signals;

[0296] If the signal within the distance clustering cluster is a translocation mutation or a complex mutation signal, directly execute steps three and four;

[0297] The edit distance calculation calls the Align function in the edlib library.

[0298] For an insertion mutation signal, the insertion mutation sequence is directly obtained from the sequence information record of the mutation signal (sequence information identifier Seq of the insertion signal);

[0299] For deletion mutation signals, the sequence with the same chromosome name as the deletion mutation signal and between the mutation start position and mutation end position of the mutation signal is extracted from the reference genome as the deletion mutation sequence;

[0300] For repeated variation signals, the sequence with the same chromosome name as the repeated variation signal and between the variation start position and variation end position of the variation signal is extracted from the reference genome as the repeated variation sequence;

[0301] For inversion mutation signals, the sequence with the same chromosome name as the inversion mutation signal and between the mutation start position and the mutation end position of the mutation signal is extracted from the reference genome as the inversion mutation sequence;

[0302] Step 3

[0303] According to the edit distance between any two signals in the obtained cluster, the following operations are performed on the signals in the cluster:

[0304] Set the similarity threshold (default 10);

[0305] If the sequence editing similarity of the inserted variant sequences of the two variant signals is less than the similarity threshold, the two variant signals are classified as one cluster; otherwise, the two variant signals are not one cluster; until all variant signals in the inserted cluster are judged;

[0306] If the sequence editing similarity of the deleted variant sequences of two variant signals is less than the similarity threshold, the two variant signals are classified as one cluster; otherwise, the two variant signals are not one cluster; until all variant signals in the deleted cluster are judged;

[0307] If the sequence editing similarity of the repeated variant sequences of two variant signals is less than the similarity threshold, the two variant signals are classified as one cluster; otherwise, the two variant signals are not one cluster; until all variant signals in the repeated cluster are judged;

[0308] If the sequence editing similarity of the inversion mutation sequences of two mutation signals is less than the similarity threshold, the two mutation signals are classified as one cluster; otherwise, the two mutation signals are not one cluster; until all the mutation signals in the inversion cluster are judged;

[0309] Execute steps three and four;

[0310] Step 3 and 4: Genotyping the variant signals within the clusters; the specific process is as follows:

[0311] If the number of variant signals in the cluster is 2 and the variant signals come from different haplotypes (the input is two bam files, so it can be determined which file the detected variant comes from, that is, which haplotype), it is determined to be the same haplotype sequence, and the genotype is 1 / 1;

[0312] If the number of variant signals in the cluster is 1, it is judged as a heterohaplotype sequence, and the genotype is judged as 0 / 1 or 1 / 0 according to the haplotype to which it belongs;

[0313] Otherwise, it is considered as noise and the corresponding cluster is discarded;

[0314] Output mutation signals containing genotyping information;

[0315] Write each mutation signal record into a file in VCF format.

[0316] There are three types of genotyping: 0 / 1, 1 / 0, and 1 / 1; 1 means that this is different from the reference genome and there is a variation; 0 means that it is the same as the reference genome and there is no variation; if it is 1 / 1, it means that the genome used for testing has two identical variation on the two haplotypes at this position; if it is 1 / 0, it means that the genome used for testing has only a variation on the first haplotype at this position; if it is 0 / 1, it means that the genome used for testing has only a variation on the second haplotype at this position.

[0317] Through the above steps, we obtain the necessary information of all mutations in this bam file (chromosome name, mutation start position, mutation end position / variation length, mutation sequence, genotype, other supporting information, etc.), and then by writing each mutation signal record into a file in VCF format, we obtain a complete VCF file with all structural variation information indicated by the bam file.

[0318] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0319] The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A variation detection method based on end-to-end assembled genome, characterized by: The specific process of the method is: Step 1: Process the variant signal to obtain the processed variant signal; the specific process is: Step 11, obtain BAM format comparison information; Step 12: Filter the BAM format comparison information obtained in step 11 to obtain filtered BAM format comparison information; Step 13: extract the variation signal from the filtered BAM format comparison information; Step 14: Cluster the extracted variant signals and output the clustering results; Step 15: Filter all the output cluster results to obtain the final cluster result, that is, integrate the signals contained in each cluster into a structural variation signal for output; Step 2: perform mixed sorting on the processed variation signals obtained in step 1, perform spatial distribution clustering on the sorted variation signals, and determine whether the variation signals after spatial distribution clustering are potential complex variations. If so, it is determined to be a complex variation signal, and a complex variation signal is obtained, and step 3 is executed; if not, it is determined to be a simple variation signal, and step 3 is executed. A simple variation signal is a variation signal of insertion, deletion, inversion, translocation, and duplication; Step 3: Perform genotyping on the variant signal obtained in step 2.

2. The method for variation detection based on end-to-end genome assembly according to claim 1, characterized in that: In the step 1, the BAM format comparison information is obtained; the specific process is: 1) Obtain sequencing data files in fastq format from a high-throughput sequencing platform; 2) Use a splicing tool to process the sequencing data file in fastq format to obtain two haplotype assembled genomes after typing. The format of the two haplotype assembled genomes is fasta format; 3) Use the alignment tool to align the two haplotype assembled genomes with the reference genome, and generate two SAM format alignment information files; 4) Convert the two SAM format comparison information files into binary BAM format comparison information files respectively.

3. The method for variation detection based on end-to-end genome assembly according to claim 2, characterized in that: In the step 12, the BAM format comparison information obtained in step 11 is filtered to obtain the filtered BAM format comparison information; the specific process is: Classify the comparison information record according to the flag bit in the BAM format comparison information to determine the type of the comparison information record; The types of comparison information records can be classified as primary comparison, secondary comparison or supplementary comparison; Through the SA tag of the main alignment, the alignment information records of all fragments constituting the same contig are obtained according to the same contig name, as given in Formula 1: contig={seg1,seg2,...,seg j ,...,say n },1<i≤n (1) The filtering conditions are: Among them, contig is the alignment information of the spliced ​​fragment, seg i It is the comparison information record of the ith fragment in the contig after the fragments are sorted according to their relative starting positions. segi Indicates seg i The end position in the contig; Indicates seg i-1 The end position in the contig; 1<i≤n, n represents the total number of fragments in the contig; Select The comparison information is used as the filtered BAM format comparison information.

4. The method for variation detection based on end-to-end genome assembly according to claim 3, characterized in that: In the step 1-3, the variation signal in the filtered BAM format comparison information is extracted; the specific process is: The comparison information is divided into intra signal and inter signal according to the flag bit in the filtered BAM format comparison information; Perform signal processing on the intra signal to obtain two types of variation signals: insertion and deletion; Perform signal processing on the inter signal to obtain five types of variation signals: insertion, deletion, inversion, translocation, and duplication; The intra signal is processed to obtain two kinds of variation signals, insertion and deletion. The specific process is as follows: If the filtered seg i If the number before the I character in the CIGAR character position in is greater than 30bp, it is considered that there is an insertion signal, and the chromosome name where the insertion signal is located, the starting position of the insertion signal on the reference genome, the insertion signal size, and the insertion signal sequence are recorded as an insertion signal; otherwise, there is no insertion signal; If the filtered seg i If the number before the D character in the CIGAR character position in is greater than 30bp, it is considered that there is a deletion signal, and the chromosome name where the deletion signal is located, the starting position of the deletion signal on the reference genome, the ending position of the deletion signal on the reference genome, and the size of the deletion signal are recorded as a deletion signal; otherwise, there is no deletion signal.

5. The method for variation detection based on end-to-end genome assembly according to claim 4, characterized in that: The inter signal is processed to obtain five types of variation signals: insertion, deletion, inversion, translocation, and duplication; the specific process is as follows: Traverse the two adjacent segs that make up the contig, denoted as seg1 and seg2 respectively, where the fragment with the smaller starting coordinate is designated as the first aligned fragment; OLP ref =RefEnd seg1 -RefStart seg2 in, Indicates the starting position of seg1 on the spliced ​​fragment, Indicates the starting position of seg2 on the spliced ​​fragment; Indicates the termination position of seg1 on the spliced ​​fragment, Indicates the termination position of seg2 on the spliced ​​fragment; Indicates the starting position of seg1 on the reference genome, Indicates the starting position of seg2 on the reference genome; Indicates the termination position of seg1 on the reference genome, Indicates the termination position of seg2 on the reference genome; read dis represents the distance between the starting position of seg2 on the spliced ​​fragment and the ending position of seg1 on the spliced ​​fragment, ref dis The OLP represents the distance between the starting position of seg2 and the ending position of seg1 on the reference genome. ref Indicates the overlap size between the end position of seg1 on the reference genome and the start position of seg2 on the reference genome; 1) For the inserted mutation signal, if the inserted mutation signal meets the following preconditions: in, Indicates the chromosome name in seg1, represents the chromosome name in seg2; Formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2; And satisfy: Then extract the insertion signal: Among them, Chrom represents the chromosome name identifier where the extracted signal is located. INS Indicates the name of the chromosome where the extracted insertion signal is located, Start indicates the starting position of the extracted signal, RefEnd seg1 Indicates the end position of the reference genome in seg1, RefEnd seg2 Indicates the end position of the reference genome in seg2, SVlen indicates the variant length, Seq indicates the sequence information identifier of the insertion signal, Seq INS Indicates the sequence corresponding to the insertion signal; If formulas (4) and (5) are not satisfied, the insertion signal is not extracted; 2) For deletion of mutation signals, if the deletion of mutation signals meets the following preconditions: in, Indicates the chromosome name in seg1, represents the chromosome name in seg2; Formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2; And satisfy: Then extract the delete signal: Among them, Chrome DEL Indicates the name of the chromosome where the extracted deletion signal is located, End indicates the end position mark of the extracted signal, and other marks are defined as above; If formulas (4) and (7) are not satisfied, the deletion signal is not extracted; 3) For repeated variant signals, if the repeated variant signals meet the following preconditions: in, Indicates the chromosome name in seg1, represents the chromosome name in seg2; Formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2; And satisfy: Then extract the repeated signal: Among them, ReadEnd seg1 Indicates the end position of the spliced ​​fragment in seg1, ReadStart seg2 Indicates the starting position of the spliced ​​fragment in seg2, RefStrart seg2 Indicates the starting position of the reference genome in seg2, Chrom DUP Indicates the name of the chromosome where the extracted DUP signal is located, Th olp Indicates the distance threshold for extracting repeated DUP signals; If formulas (4) and (9) are not satisfied, the repeated signal will not be extracted; 4) For inversion mutation signals; 41) If the inversion mutation signal meets the following preconditions: in, Indicates the chromosome name in seg1, represents the chromosome name in seg2; Formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2; And satisfy: The alignment directions in seg1 and seg2 records are different. There is a reverse alignment and a forward alignment between seg1 and seg2. The seg in the forward alignment is recorded as fwd, and the seg in the reverse alignment is recorded as rev. |ReadEnd rev -ReadStart fwd |≤100bp (11) Then extract the inversion (INV) signal: If formulas (4) and (11) are not satisfied, the inversion (INV) signal is not extracted; 42) If the inversion mutation signal meets the following preconditions: in, Indicates the chromosome name in seg1, represents the chromosome name in seg2; Formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2; And satisfy: The alignment directions in seg1 and seg2 records are different. There is a reverse alignment and a forward alignment between seg1 and seg2. The seg in the forward alignment is recorded as fwd, and the seg in the reverse alignment is recorded as rev. |ReadEnd fwd -ReadStart rev |≤100bp (13) Extract the INV signal: If formulas (4) and (13) are not satisfied, the inversion (INV) signal is not extracted; Among them, ReadEnd rev ReadStart is the end position of the spliced ​​fragment in the reverse alignment information. fwd RefStart is the starting position of the spliced ​​fragment in the forward alignment information. rev RefStart is the starting position of the reference genome in the reverse alignment information. fwd is the starting position of the reference genome in the forward alignment information, ReadEnd fwd ReadStart is the end position of the spliced ​​fragment in the forward alignment information. rev RefEnd is the starting position of the spliced ​​fragment in the reverse alignment information, fwd RefEnd is the reference genome end position in the forward alignment information. rev The reference genome end position in the reverse alignment information; Chrom INV Start indicates the chromosome name where the INV mutation signal occurs; INV Indicates the start position of the INV mutation signal, End INV The mark indicating the end position of the INV variation signal; 5) For translocation mutation signals; 51) If the translocation mutation signal meets the following preconditions: Among them, formula (15) shows that the chromosome name recorded in seg1 is inconsistent with the chromosome name recorded in seg2; And satisfy: |ReadEnd seg1 -ReadStart seg2 |≤1000bp (16) Then extract the translocation variation signal: If formulas (15) and (16) are not satisfied, the translocation signal is not extracted; 52) If the translocation mutation signal meets the following preconditions: Among them, formula (15) shows that the chromosome name recorded in seg1 is inconsistent with the chromosome name recorded in seg2; And satisfy: |ReadEnd seg2 -ReadStart seg1 |≤1000bp (18) Then extract the translocation signal: If formulas (15) and (18) are not satisfied, the translocation signal is not extracted; Among them, Chrom1 represents the identifier of the starting chromosome name where the TRA mutation signal occurs, Pos1 represents the identifier of the starting position where the TRA mutation signal occurs, Chrom2 represents the identifier of the ending chromosome name where the TRA mutation signal occurs, and Pos2 represents the identifier of the ending position where the TRA mutation signal occurs.

6. The method for variation detection based on end-to-end genome assembly according to claim 5, characterized in that: In the step 1-4, the extracted variant signals are clustered and the clustering results are output; the specific process is: (1) Targeting the extracted deletion signals, insertion signals, duplication signals, inversion signals, and translocation signals; Initialize the deletion cluster container; initialize the insertion cluster container; initialize the duplicate cluster container; initialize the inversion cluster container; initialize the translocation cluster container; (2) Sort the variant signals in each variant type set in ascending order according to the starting position where they appear in the reference genome; (3) For the signals in the mutation signal set of the deleted mutation type, perform the following operations in sequence: a. Initialize an empty cluster and add the first signal in the mutation signal set of the deletion mutation type as the starting signal, and execute d; b. Calculate the similarity between the current signal and the last signal in the current cluster; c. If the similarity is less than a predetermined threshold, the current cluster is added to the deleted cluster container and the current signal is added to the new cluster; If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster; d. Move to the next signal and repeat bc until all signals in the set of variant signals with the variant type to be deleted are determined; (4) For the signals in the mutation signal set of the inserted mutation type, perform the following operations in sequence: a. Initialize an empty cluster and add the first signal in the mutation signal set of the inserted mutation type as the starting signal, and execute d; b. Calculate the similarity between the current signal and the last signal of the current cluster; c. If the similarity is less than a predetermined threshold, the current cluster is added to the inserted cluster container and the current signal is added to the new cluster; If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster; d. Move to the next signal and repeat bc until all signals in the set of mutant signals of the inserted mutation type are determined; (5) For the signals in the mutation signal set of the repeated mutation type, perform the following operations in sequence: a. Initialize an empty cluster and add the first signal in the set of mutation signals of the repeated mutation type as the starting signal, and execute d; b. Calculate the similarity between the current signal and the last signal of the current cluster; c. If the similarity is less than a predetermined threshold, the current cluster is added to the repeated cluster container and the current signal is added to the new cluster; If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster; d. Move to the next signal and repeat bc until all signals in the set of variant signals of the repeated variant type are determined; (6) For the signals in the mutation signal set of the inversion mutation type, the following operations are performed in sequence: a. Initialize an empty cluster, add the first signal in the mutation signal set of the inversion mutation type as the starting signal, and execute d; b. Calculate the similarity between the current signal and the last signal of the current cluster; c. If the similarity is less than a predetermined threshold, the current cluster is added to the inverted cluster container, and the current signal is added to the new cluster; If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster; d. Move to the next signal and repeat bc until all signals in the set of mutation signals of the inversion mutation type are determined; (7) For the signals in the mutation signal set of the translocation mutation type, the following operations are performed in sequence: a. Initialize an empty cluster, add the first signal in the set of mutation signals of the translocation mutation type as the starting signal, and execute d; b. Calculate the similarity between the current signal and the last signal of the current cluster; c. If the similarity is less than a predetermined threshold, the current cluster is added to the translocation cluster container, and the current signal is added to the new cluster; If the similarity is greater than or equal to a predetermined threshold, the current signal is added to the current cluster; d. Move to the next signal and repeat bc until all signals in the set of variant signals of the translocation variant type are determined; (8) Obtain all clustering results of all variant type signals.

7. The method for variation detection based on end-to-end genome assembly according to claim 6, characterized in that: In b, the similarity between each signal and the last signal of the current cluster is calculated; the specific process is:

1. For insertion, deletion, inversion, and duplication mutation signals: 1) Calculate the location similarity S pos Scoring; the specific process is: Position similarity S pos The score calculation formula is as follows: S pos =|Start1-Start2| (20) Among them, Start1 and Start2 are the starting positions of the signals in the two variant signals respectively; 2) Calculate the similarity score S of the mutation signal size size ; The specific process is: The mutation lengths of the two mutation signals are SVlen1 and SVlen2 respectively, and the similarity score of the mutation length of the mutation signal is S size Calculated by the following formula: 3) Calculate the chromosome similarity score S Chr ; The specific process is: Assume that the chromosome numbers of the two mutation signals are Chrom1 and Chrom2, and the chromosome similarity score is S Chr Calculated by the following discrete function: When the chromosome names of the two variant signals are the same, the chromosome similarity score is 1; When the chromosome names of the two variant signals are different, the chromosome similarity score is 0; 4) The position similarity S pos Score, variation size similarity score S size Score, chromosome similarity score S Chr Perform weighted summation to obtain a comprehensive similarity score S; expressed as: S=S Chr ×(w size ×S size +w pos ×S pos ) (23) Among them, w pos 、w size are the weight coefficients of the position similarity score and the variation size similarity score, respectively, and w pos +w size =1; 2. For translocation mutation signals: Among them, S Chr1 represents the starting chromosome similarity score, Indicates the name of the starting chromosome of the first signal, indicates the name of the starting chromosome of the second signal; Among them, S Chr2 represents the termination chromosome similarity score, The name of the chromosome that terminates the first signal, The name of the chromosome that indicates the termination of the second signal; The starting chromosome similarity score S Chr1 , end chromosome similarity score S Chr2 Perform weighted summation to obtain a comprehensive similarity score S; expressed as: in, refers to the starting position of the second signal, Refers to the starting position of the first signal, Refers to the end position of the second signal Refers to the end position of the first signal; w pos1 is the weight coefficient of the starting position, w pos2 is the weight coefficient of the end position.

8. The method for variation detection based on end-to-end genome assembly according to claim 7, characterized in that: In the step 15, all the output cluster results are filtered to obtain the final cluster result, that is, the signals contained in each cluster are integrated into a structural variation signal for output; the specific process is: 1) Set the maximum number of clusters max_cluster; 2) If the number of clusters in a cluster is greater than max_cluster, the corresponding cluster is discarded; 3) 31) Clustering results for insertion mutation types: For each clustering cluster of the clustering result of the inserted mutation type, if the number of clusters in the clustering cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and average length of all mutation data in the clustering cluster; the specific process is: The mutation start positions of each mutation signal in the δth cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ = 1, 2, ..., N, N represents the total number of clusters in the clustering results of the insertion mutation type; The lengths of all variants in the δth cluster are summed up and divided by the number of variants to calculate the average length; Keep the other information in the δth cluster the same as the first mutation signal in the cluster; The other information is the chromosome name, the name from the contig, and the inserted sequence; The average mutation start position, average length and other information are output as the δth clustering result; Until the output results of all clustering clusters of the clustering results of the insertion mutation type are obtained, that is, the final insertion mutation signal; 32) Clustering results for deletion mutation types: For each clustering cluster of the clustering results of the deleted variant type, if the number of clusters in the clustering cluster is less than or equal to max_cluster, calculate the average variant signal starting position, average variant signal ending position and average length of all variant data in the clustering cluster; the specific process is: The mutation start positions of each mutation signal in the δ′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ′=1,2,…,N′, N′ represents the total number of clusters in the clustering result of deleting the mutation type; The mutation termination positions of each mutation signal in the δ′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position; The lengths of all variants in the δ′th cluster are summed up and divided by the number of variants to calculate the average length. Keep the other information in the δ′th cluster the same as the first mutation signal in the cluster; The other information is the chromosome name and the name from the contig; The average mutation start position, average mutation end position, average length and other information are output as the δ′th cluster; Until the output results of all clustering clusters of the clustering results of the deleted mutation type are obtained, that is, the final deleted mutation signal; 33) Clustering results for inversion mutation types: For each cluster of the inversion mutation type clustering result, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and the average mutation signal ending position of all mutation data in the cluster; the specific process is: The mutation start positions of each mutation signal in the δ″th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ″=1,2,…,N″, N″ represents the total number of clusters in the clustering results of the inversion mutation type; The mutation termination positions of each mutation signal in the δ″th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position; Keep the other information in the δ″th cluster the same as the first mutation signal in the cluster; The other information is the chromosome name and the name from the contig; The average mutation start position, average mutation end position and other information are used as the output result of the δ″th cluster; Until the output results of all clustering clusters of the clustering results of the inversion mutation type are obtained, that is, the final inversion mutation signal; 34) Clustering results for translocation variant types: For each cluster of the clustering results of the translocation variant type, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and the average mutation signal ending position of all mutation data in the cluster; the specific process is: The mutation start positions of each mutation signal in the δ″′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; δ″′=1,2,…,N″′, N″′ represents the total number of clusters in the clustering results of the translocation mutation type; The mutation termination positions of each mutation signal in the δ″′th cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position; Keep the other information in the δ″′th cluster the same as the first mutation signal in the cluster; The other information is the name of the starting chromosome, the name of the ending chromosome, and the name from the contig; The average mutation start position, average mutation end position and other information are output as the δ″′th cluster; Until the output results of all clustering clusters of the clustering results of the translocation variation type are obtained, that is, the final translocation variation signal; 35) Clustering results for repeated variation types: For each cluster of the clustering results of the repeated mutation type, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average mutation signal starting position and the average mutation signal ending position of all mutation data in the cluster; the specific process is: The first The mutation start positions of each mutation signal in each cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation start position; N″′ represents the total number of clusters in the clustering results of repeated variation types; The first The mutation termination positions of each mutation signal in each cluster are accumulated one by one, and then divided by the number of mutations to obtain the average mutation termination position; The first The other information in the cluster remains the same as the first mutation signal in the cluster; The other information is the chromosome name and the name from the contig; The average mutation start position, average mutation end position and other information are used as the The cluster output results are: Until the output results of all clustering clusters of the clustering results of the repeated variation type are obtained, that is, the final repeated variation signal.

9. The method for variation detection based on end-to-end genome assembly according to claim 8, characterized in that: In the step 2, the processed variation signals obtained in the step 1 are mixed and sorted, and the sorted variation signals are spatially distributed and clustered to determine whether the variation signals after spatial distribution and clustering are potential complex variations. If so, it is determined to be a complex variation signal, and a complex variation signal is obtained, and step 3 is executed; if not, it is determined to be a simple variation signal, and step 3 is executed. The simple variation signal is a variation signal of insertion, deletion, inversion, translocation, and duplication. The specific process is: 1) Use merge sort or quick sort to sort the processed variant signals obtained in step 1 in ascending order according to their starting positions on the reference genome; 2) Taking the genome position where the first mutation signal after sorting is located as the center, extend Dist / 2 base pairs forward from the center, and extend Dist / 2 base pairs backward from the center to form a Dist scanning window of length; Acquire all mutation signals within the scanning window; 3) Determine the mutation type to which all mutation signals in the scanning window belong; The types of variants include insertion, deletion, inversion, translocation, or duplication; When the number of variant signal types in a certain scanning window exceeds two, the scanning window is identified as a potential complex variant region; otherwise, the scanning window is not a potential complex variant region; 4) Extract detailed information of each complex variation signal in the potential complex variation region. The detailed information includes the variation type, variation chromosome name, and variation start position of each variation signal. The format of the complex variation signal is (CSV, chromosome name, variation type, variation start position); Among them, CSV is the format identifier; 5) Move the scanning window back one variant signal in order, and repeat 1)-4) until all complex variant signals after mixed sorting are traversed and processed.

10. The method for variation detection based on end-to-end genome assembly according to claim 9, characterized in that: In step three, genotyping is performed on the variation signal obtained in step two; the specific process is as follows: Step 3.1: Distance clustering is performed on different types of variant signals obtained in step 2. The specific process is as follows: For insertion type mutation signals, three conditions are checked: first, the chromosome name where the current mutation signal is located is the same as the chromosome name where the mutation is located in the current cluster; second, the Euclidean distance between the chromosome starting position in the current mutation signal and the chromosome starting position of the last mutation signal in the current cluster is less than 200bp; third, the difference between the mutation length of the current mutation signal and the mutation length of the last mutation signal in the current cluster is less than 50bp; if all three conditions are met, the current mutation signal is added to the current cluster; otherwise, a new cluster is opened and the current mutation signal is added to the new cluster; For deletion, inversion or duplication type mutation signals, three conditions are checked: first, the chromosome name where the current mutation signal is located is the same as the chromosome name where the mutation is located in the current cluster; second, the Euclidean distance between the mutation start position of the current mutation signal and the mutation start position of the last mutation signal in the current cluster is less than 200bp; third, the Euclidean distance between the mutation end position in the current mutation signal and the mutation end position in the last mutation signal in the current cluster is less than 200bp; if all three conditions are met, the current mutation signal is added to the current cluster; otherwise, a new cluster is opened and the current mutation signal is added to the new cluster; For the mutation signal of the translocation type, four conditions are checked: first, the chromosome name where the current mutation signal starts is the same as the chromosome name where the last mutation signal starts in the current cluster; second, the chromosome name where the current mutation signal ends is the same as the chromosome name where the last mutation signal ends in the current cluster; third, the Euclidean distance between the mutation start position of the current mutation signal and the mutation start position of the last mutation signal in the current cluster is less than 200bp; fourth, the Euclidean distance between the mutation end position of the current mutation signal and the mutation end position of the last mutation signal in the current cluster is less than 200bp; if all four conditions are met, the current mutation signal is added to the current cluster; otherwise, a new cluster is opened and the current mutation signal is added to the new cluster; For the mutation signal of complex mutation, three conditions are checked: first, the chromosome name where the current mutation signal starts is the same as the chromosome name where the last mutation starts in the current cluster; second, the Euclidean distance between the mutation start position of the current mutation signal and the last mutation start position in the current cluster is less than 1000bp; third, the mutation type of the current mutation signal is the same as the mutation type of the last mutation signal in the current cluster; if all three conditions are met, the current signal is added to the current cluster; otherwise, a new cluster is opened and the current signal is added to the new cluster; Step 32 If the signal within the distance clustering cluster is an insertion mutation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the insertion mutation sequences of any two mutation signals; If the signal within the distance clustering cluster is a deletion mutation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the deletion mutation sequences of any two mutation signals; If the signal within the distance clustering cluster is a repeated variation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the repeated variation sequences of any two variation signals; If the signal within the distance clustering cluster is an inversion mutation signal, the minimum edit distance algorithm is used to calculate the sequence editing similarity between the inversion mutation sequences of any two mutation signals; If the signal within the distance clustering cluster is a translocation mutation or a complex mutation signal, directly execute steps three and four; Step 3 Set similarity threshold; If the sequence editing similarity of the inserted variant sequences of two variant signals is less than the similarity threshold, the two variant signals are classified into one cluster; Otherwise, the two variant signals are not a cluster; until all variant signals in the inserted cluster are judged; If the sequence editing similarity of the deletion variant sequences of two variant signals is less than the similarity threshold, the two variant signals are classified into one cluster; Otherwise, the two mutation signals are not a cluster; Until all the variant signals in the cluster are deleted; If the sequence editing similarity of the repeated variant sequences of two variant signals is less than the similarity threshold, the two variant signals are classified as one cluster; otherwise, the two variant signals are not one cluster; until all variant signals in the repeated cluster are judged; If the sequence editing similarity of the inversion mutation sequences of two mutation signals is less than the similarity threshold, the two mutation signals are classified as one cluster; otherwise, the two mutation signals are not one cluster; until all the mutation signals in the inversion cluster are judged; Execute steps three and four; Step 3 and 4: Genotyping the variant signals within the clusters; the specific process is as follows: If the number of variant signals in a cluster is 2 and the variant signals come from different haplotypes, they are judged to be the same haplotype sequence and the genotype is 1 / 1; If the number of variant signals in the cluster is 1, it is judged as a heterohaplotype sequence, and the genotype is judged as 0 / 1 or 1 / 0 according to the haplotype to which it belongs; Otherwise, it is considered as noise and the corresponding cluster is discarded; Output mutation signals containing genotyping information; Write each mutation signal record into a file in VCF format.

Images