Variant detection method based on end-to-end assembly of genomes

Abstract

Based on end-to-end assembly genome-based variation detection method, the present application relates to a variation detection method. The purpose of the present application is to solve the problems of incomplete detection, inaccurate detection, inability to identify complex variation types and low accuracy of genotyping in the current method based on splicing detection structural variation. The process of the end-to-end assembly genome-based variation detection method is as follows: step one, variation signal processing, obtaining the processed variation signal; step two, mixed sorting of the processed variation signal obtained in step one, spatial distribution clustering of the sorted variation signal, judging whether the spatial distribution clustered variation signal is a potential complex variation, if yes, judging the complex variation type and outputting, obtaining the complex variation signal, and executing step three; if not, it is a simple variation signal, and step three is executed; step three, genotyping of the variation signal obtained in step two. The present application is used in the field of variation detection.

Description

TECHNICAL FIELD

[0001] The present application relates to a variation detection method. BACKGROUND

[0002] At present, the means for analyzing the gene sequence of an organism is to obtain a DNA sequence fragment with high reliability through sequencing and quality control, and then analyze the obtained fragment according to the corresponding task through various tools and methods. The most important analysis is to obtain the difference from the known reference sequence, that is, the variation in biology; the variation in a long and 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 great impact on the traits of an organism, so it is very important to detect structural variation. At present, the common process for detecting structural variation on the human genome is mostly based on second-generation or third-generation sequencing data, but due to the perfection of splicing tools and technological progress, it is possible to produce individual typing of end-to-end spliced chromosomes. Existing research also shows that structural variation detection based on splicing of genomes can find variations that cannot be found by second-generation and third-generation sequencing, and has an advantage in detecting complex structural variations. However, the current method for detecting structural variation based on splicing has the problems of incomplete detection, inaccurate detection, inability to identify complex variation types, and low accuracy of gene typing.

[0003] Therefore, the present method provides a process and method for detecting structural variation based on end-to-end splicing, which improves the accuracy, sensitivity and detection efficiency of existing genome structural variation detection. SUMMARY

[0004] The purpose of the present application is to solve the problems of incomplete detection, inaccurate detection, inability to identify complex variation types and low accuracy of gene typing in the current method for detecting structural variation based on splicing, and to propose a variation detection method based on end-to-end assembled genomes.

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

[0006] Step one, variation signal processing to obtain processed variation signals; the specific process is as follows:

[0007] Step one, obtaining BAM format alignment information;

[0008] Step two, filtering the BAM format alignment information obtained in step one to obtain filtered BAM format alignment information;

[0009] Step three, extracting variation signals from the filtered BAM format alignment information;

[0010] Step four, clustering the extracted variant signals, and outputting the clustering cluster results;

[0011] Step five, filtering all the output clustering cluster results, and obtaining the final clustering results, that is, the signals contained in each clustering cluster are integrated into a structural variant signal for output;

[0012] Step two, mixing and sorting the processed variant signals obtained in step one, and performing spatial distribution clustering on the sorted variant signals; judging whether the variant signals after spatial distribution clustering are potential complex variants; if yes, judging that they are complex variant signals, and obtaining the complex variant signals, and executing step three; if not, judging that they are simple variant signals, and executing step three; the simple variant signals are one of insertion, deletion, inversion, translocation and duplication;

[0013] Step three, genotyping the variant signals obtained in step two.

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

[0015] The application provides a set of end-to-end assembled genome-based variation detection method flow and method, which divides the variation detection flow into three parts: a variation signal processing part, a complex structure variation detection part and a gene typing part, and can realize efficient and reliable detection of the structure variation of the spliced genome. In the variation signal processing, the alignment information filtering method of maximizing non-overlapping alignment information is screened out to reduce false positives in the variation detection; in the variation signal extraction, a rule scheme for detecting basic structure variation is given, different threshold parameters are set for different variation types, so that the variation signal is more accurately and comprehensively extracted; in the variation signal clustering and filtering, the unsupervised clustering and filtering method based on feature similarity is proposed, the similarity function is used to measure the similarity of the variation signal, so as to realize the clustering of similar signals in each variation type and form the average variation to achieve the filtering effect; at the same time, the abnormal variation cluster signal is discarded, and the credibility of the detected variation signal is improved. In the complex structure variation detection, the potential complex variation site is identified, and then the complex structure variation fusion is carried out to realize the detection of the complex structure variation type. In the gene typing part, the gene typing method based on distance similarity clustering and sequence similarity classification is proposed, the distance clustering and hierarchical clustering are combined, and the minimum edit distance is used to more finely divide the haplotype, so as to solve the precision and efficiency problem of the structure variation typing. The application can accurately and efficiently detect common structure variations and complex structure variation types of the genome and accurately and efficiently complete the gene typing function by providing a complete end-to-end spliced genome structure variation detection flow and method, and fills the gap and deficiency of the existing method in the end-to-end genome structure variation detection. At the same time, the clustering and filtering method proposed in the application can effectively reduce false positives, avoid the influence of noise signals and improve the detection accuracy. The complex variation structure variation detection method can quickly identify the complex structure variation and provide clues for revealing the manifestation form and rule of the complex variation mode. The gene typing method effectively solves the precision and efficiency problem of the structure variation typing. The overall scheme is fast, efficient and accurate, and has good application value. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 is the overall flowchart of the application;

[0017] Figure 2 is the variation signal clustering flowchart;

[0018] Figure 3 is the gene typing flowchart. DETAILED DESCRIPTION

[0019] Specific implementation one: the specific process of the embodiment based on the end-to-end assembled genome variation detection method is as follows:

[0020] To get the final variant file from the assembled genome alignment file, the algorithm mainly contains three steps: variant signal processing part, complex structure variant detection part and genotyping part.

[0021] To obtain the bam file for analysis, first obtain the quality-controlled sequencing data fastq format file from high-throughput sequencing platforms such as Illumina, PacBio, etc., then assemble and genotype two haplotype genome assembly fasta files by assembly tools such as Hifiasm, Peregrine, Wtdbg2, etc., then align the assembled genome with the reference genome using alignment tools such as minimap2, NGMLR, Winnowmap and LRA to generate SAM format files, and then convert SAM to binary BAM files with samtools.

[0022] In the variant signal processing, the overall process is to get all the alignment information records from the bam file, and the information in the alignment record is divided into intra signal and inter signal according to its flag bit, and the signal processing is carried out respectively to obtain five kinds of basic variant signals.

[0023] To realize the processing of variant signals, it is specifically divided into three steps: alignment information filtering, variant signal extraction and variant signal clustering.

[0024] Step one, variant signal processing, get the processed variant signal; the specific process is:

[0025] Step one, get the BAM format alignment information;

[0026] Step two, filter the BAM format alignment information obtained in step one to obtain the filtered BAM format alignment information;

[0027] Step three, extract the variant signal in the filtered BAM format alignment information;

[0028] Step four, cluster the extracted variant signal, and output the clustering cluster result;

[0029] Step five, filter all the output clustering cluster results to get the final clustering result, that is, integrate the signals contained in each clustering cluster into a structure variant signal for output;

[0030] Step two, mix and sort the processed variant signals obtained in step one, cluster the sorted variant signals in space distribution, and determine whether the variant signals clustered in space distribution are potential complex variants. If yes, it is determined as a complex variant signal, and a complex variant signal is obtained, and step three is executed. If not, it is determined as a simple variant signal, and step three is executed. The simple variant signal is one of insertion, deletion, inversion, translocation and duplication;

[0031] Step three, genotype the variant signals obtained in step two.

[0032] Specific implementation method two: the difference between this embodiment and specific implementation method one is that the BAM format alignment information is obtained in step one. The specific process is as follows:

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

[0034] Obtain fastq format sequencing data files from an Illumina, PacBio high-throughput sequencing platform;

[0035] 2) Process the fastq format sequencing data files by using a splicing tool to obtain two haplotype assembled genomes after genotyping, and the format of the two haplotype assembled genomes is fasta format;

[0036] Genotyping refers to that the genome of a diploid organism is usually one from the father and one from the mother, and it is necessary to determine whether the gene is from the father or the mother;

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

[0038] 3) Align the two haplotype assembled genomes with the reference genome respectively by using an alignment tool to generate two SAM format alignment information files;

[0039] Align the assembled genome with the reference genome by using an alignment tool such as minimap2, NGMLR, Winnowmap and LRA to generate a SAM format file;

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

[0041] Convert the SAM format alignment information into binary BAM format alignment information by using samtools.

[0042] The other steps and parameters are the same as those in specific implementation method one.

[0043] Specific implementation three: the difference between this embodiment and specific implementation one or two is that the BAM format alignment information obtained in step one is filtered in step two to obtain filtered BAM format alignment information; the specific process is as follows:

[0044] According to the flag flag in the BAM format alignment information, the alignment information records are classified to determine the type of the alignment information records;

[0045] The BAM format alignment information is composed of many lines, and each line represents a record;

[0046] The type of the alignment information record can be primary alignment, secondary alignment, or supplementary alignment, etc.

[0047] The read / contig name of each alignment information record itself, the chromosome name of the alignment, the start position of the alignment to the reference genome, the alignment quality, the CIGAR string of the alignment, the sequence of the alignment, the supplementary information of the alignment, etc.

[0048] The supplementary information contains other tags to provide additional information, such as the SA tag recording all split alignment information constituting a same contig fragment, that is, according to the same contig name (the first column), all alignment information records of all fragments constituting the same contig are obtained through the SA tag of the primary alignment, as shown in formula 1:

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

[0050] The fasta file is input into the alignment software to obtain the bam file, and the alignment software can split the same name fasta record into several alignment records, and the several alignment records come from the same fasta record, that is, from the same contig splicing fragment; the alignment refers to finding the position of the gene sequence in the reference genome, in order to clearly indicate the position, it is necessary to record which chromosome (chromosome name) the gene sequence is aligned to, and which position (chromosome start position) on the chromosome the gene sequence is aligned to, in order to indicate the credibility of the alignment, there is an alignment quality to measure, the higher the quality, the better and more reliable the alignment; the CIGAR string of the alignment is to indicate the form of the alignment, and the alignment sequence is restored to the reference genome sequence through what operation, more (insertion, I), less (deletion, D), or consistent (M) operation; the alignment sequence refers to the specific sequence AGCT of the alignment; the gene sequence is formed by the arrangement and combination of four characters AGCT;

[0051] Before the conventional variant signal processing, there is usually a simple comparison of the alignment quality of the alignment information, and the alignment information with alignment quality greater than a certain threshold (generally 20) is screened out; but since there is a redundant alignment phenomenon in the alignment information, that is, the same gene sequence may be aligned to multiple positions of the genome to generate multiple alignment information records, therefore, a filtering condition between alignment information records is provided here to obtain maximized non-overlapping alignment fragments to reduce false positives;

[0052] The filtering condition is:

[0053]

[0054] Wherein, contig is the alignment information of the spliced fragment, seg i is the alignment information record of the i-th segment in contig sorted according to the relative position of the start of the segment, ReadEnd segi represents the end position of seg i in contig; represents the end position of seg i-1 in contig; 1 < i < n, n represents the total number of segments in contig;

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

[0056] The other steps and parameters are the same as those in embodiment one or two.

[0057] Embodiment four: different from one of embodiments one to three, the step one to three extracts the variant signal from the filtered BAM format alignment information; the specific process is:

[0058] According to the flag flag in the filtered BAM format alignment information, the alignment information is divided into intra signal and inter signal;

[0059] The intra signal is processed to obtain two kinds of variant signals of insertion and deletion;

[0060] The inter signal is processed to obtain five kinds of variant signals of insertion, deletion, inversion, translocation and duplication;

[0061] The intra signal is processed to obtain two kinds of variant signals of insertion and deletion; the specific process is:

[0062] If the filtered seg iIf the number before the I character in the CIGAR character position in the is greater than 30bp, it is considered that there is an insertion signal, and the chromosome name, 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 the is greater than 30bp, it is considered that there is a deletion signal, and the chromosome name, 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 insertion and deletion signals greater than 30bp directly obtained from the CIGAR character position in the alignment information record of each filtered seg i The format of CIGAR is a number plus a 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 the chromosome name, the starting position on the reference genome, the insertion size, and the insertion sequence are recorded as an insertion signal. For example, “60I” indicates that there is an insertion signal with a length of 60, and the remaining 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 the chromosome name, the starting position on the reference genome, the ending position on the reference genome, and the deletion size are recorded as a deletion signal.

[0065] The other steps and parameters are the same as one of the first to third embodiments.

[0066] Embodiment five: The difference between this embodiment and one of the first to fourth embodiments is that the inter signal is subjected to signal processing to obtain five kinds of variation signals of insertion, deletion, inversion, translocation, and duplication. The specific process is as follows:

[0067] The Inter signal is set by filtering the information recorded between segs about the difference in the distance before and after the contig, and the difference in the chromosome number and physical distance aligned on the reference genome, and implementing heuristic rules. For example, for the five basic types of variation signals (insertion, deletion, inversion, translocation, and duplication), different distance threshold judgment criteria are formulated according to experience and previous data analysis results: if a signal is in a specific range in the distance before and after the contig, and is aligned to a specific chromosome on the reference genome and the distance meets another condition, it is determined to meet the extraction requirements, triggering the extraction process. Specifically, five specific variation signals are extracted from read dis and ref dis and the corresponding information, as follows:

[0068] Two adjacent segs constituting the contig are traversed, denoted as seg1 and seg2, where the segment with a smaller starting coordinate is designated as the first aligned segment;

[0069]

[0070] wherein, represents the starting position of seg1 on the contig, represents the starting position of seg2 on the contig;

[0071] represents the end position of seg1 on the contig, represents the end position of seg2 on the contig;

[0072] represents the starting position of seg1 on the reference genome, represents the starting position of seg2 on the reference genome;

[0073] represents the end position of seg1 on the reference genome, represents the end position of seg2 on the reference genome;

[0074] read dis represents the distance between the starting position of seg2 on the contig and the end position of seg1 on the contig, ref dis represents the distance between the starting position of seg2 on the reference genome and the end position of seg1 on the reference genome, OLP ref represents the overlap size of the end position of seg1 on the reference genome and the starting position of seg2 on the reference genome on the reference genome;

[0075] 1) For an insertion (INS) variant signal, if the insertion (INS) variant signal satisfies the following preconditions:

[0076]

[0077] wherein, denotes the chromosome name in seg1, denotes the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2;

[0078] and satisfies:

[0079]

[0080] then the insertion (INS) signal is extracted:

[0081]

[0082] wherein, Chrom denotes the chromosome name identifier where the extracted signal is located, Chrom INS denotes 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 denotes the start position identifier of the extracted signal, RefEnd seg1 denotes the end position of the reference genome in seg1, RefEnd seg2 denotes the end position of the reference genome in seg2, SVlen denotes the length of the variant, Seg denotes the sequence information identifier of the insertion signal, Seq INS denotes the sequence corresponding to the insertion signal;

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

[0084] 2) For a deletion (DEL) variant signal, if the deletion (DEL) variant signal satisfies the following preconditions:

[0085]

[0086] wherein, denotes the chromosome name in seg1, denotes the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2;

[0087] and satisfies:

[0088]

[0089] then the deletion (DEL) signal is extracted:

[0090]

[0091] where Chrom DEL denotes the chromosome name of the extracted deletion signal (chromosome name recorded in the first alignment information or the second alignment information), End denotes the end position identifier of the extracted signal, and other identifiers are defined as above;

[0092] If formula (4), (7) is not satisfied, the deletion (DEL) signal is not extracted;

[0093] 3) For the duplication (DUP) variation signal, if the duplication (DUP) variation signal satisfies the following precondition:

[0094]

[0095] where, denotes the chromosome name in seg1, denotes the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2;

[0096] and satisfies:

[0097]

[0098] the duplication (DUP) signal is extracted:

[0099]

[0100] where, ReadEnd seg1 denotes the end position of the spliced fragment in seg1, ReadStart seg2 denotes the start position of the spliced fragment in seg2, RefStrart seg2 denotes the start position of the reference genome in seg2, Chrom DUP denotes the chromosome name of the extracted DUP signal (chromosome name recorded in the first alignment information or the second alignment information), Th olp denotes the distance threshold value for extracting the duplication DUP signal (the default value is 5000bp), and different distance threshold values are set according to different types of variations, so as to more accurately identify each variation;

[0101] If formula (4), (9) is not satisfied, the duplication (DUP) signal is not extracted;

[0102] 4) For the inversion (INV) variation signal;

[0103] 41) If the inversion (INV) variation signal satisfies the following precondition:

[0104]

[0105] wherein, represents the chromosome name in seg1, represents the chromosome name in seg2; formula (4) means that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2;

[0106] and satisfy:

[0107] the alignment directions in seg1 and seg2 records are different, seg1 and seg2 exist one reverse alignment and one forward alignment, the seg of forward alignment is recorded as fwd, and the seg of reverse alignment is recorded as rev;

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

[0109] then the inversion (INV) signal is extracted:

[0110]

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

[0112] 42), if the inversion (INV) variation signal satisfies the following precondition:

[0113]

[0114] wherein, represents the chromosome name in seg1, represents the chromosome name in seg2; formula (4) means that the chromosome name recorded in seg1 is consistent with the chromosome name recorded in seg2;

[0115] and satisfy:

[0116] the alignment directions in seg1 and seg2 records are different, seg1 and seg2 exist one reverse alignment and one forward alignment, the seg of forward alignment is recorded as fwd, and the seg of reverse alignment is recorded as rev;

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

[0118] extract the INV signal:

[0119]

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

[0121] wherein, ReadEnd rev is the end position of the splice in the forward alignment information, ReadStart fwd is the start position of the splice in the forward alignment information, RefStart rev is the start position of the reference genome in the reverse alignment information, RefStart fwd is the start position of the reference genome in the forward alignment information, ReadEnd fwd is the end position of the splice in the forward alignment information, ReadStart rev is the start position of the splice in the reverse alignment information, RefEnd fw is the end position of the reference genome in the forward alignment information, RefEnd rev is the end position of the reference genome in the reverse alignment information; Chrom INV represents the identifier of the chromosome name where the INV variation signal occurs; Start INV represents the identifier of the start position of the INV variation signal, End INV represents the identifier of the end position of the INV variation signal;

[0122] 5) for translocation (TRA) variation signal;

[0123] 51) if the translocation (TRA) variation signal satisfies the following precondition:

[0124]

[0125] wherein, formula (15) indicates that the chromosome name recorded in seg1 is inconsistent with the chromosome name recorded in seg2;

[0126] and satisfies:

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

[0128] then the TRA signal is extracted:

[0129]

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

[0131] 52) if the translocation (TRA) variation signal satisfies the following precondition:

[0132]

[0133] wherein formula (15) indicates that the chromosome name recorded in segl is inconsistent with the chromosome name recorded in seg2;

[0134] and satisfy:

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

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

[0137]

[0138] if formula (15), (18) is not satisfied, then do not extract the translocation (TRA) signal;

[0139] wherein Chroml represents the identifier of the starting chromosome name of the TRA variant signal, Posl represents the identifier of the starting position of the TRA variant signal, Chrom2 represents the identifier of the ending chromosome name of the TRA variant signal, and Pos2 represents the identifier of the ending position of the TRA variant signal.

[0140] The other steps and parameters are the same as one of the first to fourth embodiments.

[0141] Embodiment six: the difference between this embodiment and one of the first to fifth embodiments is that the extracted variant signals are clustered in steps one to four, and the clustering result is output; the specific process is as follows:

[0142] The existence of variant signal clustering is necessary. There are two reasons: first, due to the bias of the software for assembling genomes and aligning genome sequences and the existence of a large number of repetitive regions in the genome itself, there are multiple same type of variants in a certain interval, and most of them are likely to be false positives; second, there are some same type of variants that are distributed close to each other, which can be aggregated into one variant. Here, an unsupervised clustering and filtering method based on feature similarity is provided: for each type of variant, based on the distance distribution, the variants with similar features and close to each other are clustered together to form an average variant to achieve the filtering effect. For the same type of variant with an abnormal prominent aggregation site, it is selected to be discarded, because in the current technical flow, the credibility is not high, and the false positive probability is increased.

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

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

[0145] initialize a cluster container, and create a corresponding cluster container for each of the predetermined mutation types (e.g., DEL, INS, DUP, INV, and TRA) for the input set of mutation features;

[0146] (2) sort the mutation signals in each set of mutation types in ascending order according to the starting positions in the reference genome;

[0147] sort the sets of mutation signals for each mutation type according to a predetermined sorting rule to sort the mutations in ascending order according to the starting positions in the reference genome;

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

[0149] a. initialize an empty cluster, and add the first signal in the set of deletion mutation signals as a starting signal, and perform 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, add the current cluster to the deletion (DEL) cluster container, and add the current signal to a new cluster;

[0152] if the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster;

[0153] d. move to the next signal, and repeat b-c until all signals in the set of deletion (DEL) mutation signals are judged;

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

[0155] a. initialize an empty cluster, and add the first signal in the set of insertion mutation signals as a starting signal, and perform d;

[0156] b. calculate the similarity between the current signal and the last signal in the current cluster;

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

[0158] if the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster;

[0159] d. move to the next signal, repeat b-c until all signals in the set of variant signals of the INS variant type are judged;

[0160] (5) for signals in the set of variant signals of the DUP variant type, perform the following operations in turn:

[0161] a. initialize an empty cluster and add the first signal in the set of variant signals of the DUP variant type as a starting signal, and perform d;

[0162] b. calculate the similarity of the current signal and the last signal in the current cluster;

[0163] c. if the similarity is less than a predetermined threshold, add the current cluster to the DUP cluster container and add the current signal to a new cluster;

[0164] if the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster;

[0165] d. move to the next signal, repeat b-c until all signals in the set of variant signals of the DUP variant type are judged;

[0166] (6) for signals in the set of variant signals of the INV variant type, perform the following operations in turn:

[0167] a. initialize an empty cluster and add the first signal in the set of variant signals of the INV variant type as a starting signal, and perform d;

[0168] b. calculate the similarity of the current signal and the last signal in the current cluster;

[0169] c. if the similarity is less than a predetermined threshold, add the current cluster to the INV cluster container and add the current signal to a new cluster;

[0170] if the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster;

[0171] d. move to the next signal, repeat b-c until all signals in the set of variant signals of the INV variant type are judged;

[0172] (7) for signals in the set of variant signals of the TRA variant type, perform the following operations in turn:

[0173] a. initialize an empty cluster and add the first signal in the set of variant signals of the TRA variant type as a starting signal, and perform 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, add the current cluster to the translocation (TRA) cluster container and add the current signal to the new cluster;

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

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

[0178] (8) Obtain all cluster results of all signals of all variant types.

[0179] The other steps and parameters are the same as one of embodiments one to five.

[0180] Embodiment seven: The embodiment is different from one of embodiments one to six in that the similarity between each signal and the last signal of the current cluster is calculated in b. The specific process is as follows:

[0181] I. For insertion, deletion, inversion, and duplication variant signals:

[0182] 1) Calculate the position similarity S pos The specific process is as follows:

[0183] The position similarity S pos is calculated by the following formula:

[0184] S pos = |Start1-Start2| (20)

[0185] Wherein, Start1 and Start2 are the starting positions of the two variant signals, and the position difference is extremely important in similarity evaluation.

[0186] 2) Calculate the size similarity score S size of the variant signal; the specific process is as follows:

[0187] Let the variant length of the two variant signals be SVlen1 and SVlen2 respectively, and the size similarity score S size of the variant length of the variant signal is calculated by the following formula:

[0188]

[0189] For inversion and duplication variant signals, SVlen can be calculated by the difference between the variant termination position and the variant starting position.

[0190] This simple ratio can intuitively reflect the degree of fit of two mutation sizes, the value range is between 0 and 1, the closer to 1, the more similar the two sizes.

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

[0192] Suppose the chromosome numbers of the two mutation signals are Chrom1 and Chrom2 respectively, and the chromosome similarity score S Chr is calculated by the following discrete function:

[0193]

[0194] The chromosome similarity score is 1 when the chromosome names of the two mutation signals are the same;

[0195] The chromosome similarity score is 0 when the chromosome names of the two mutation signals are different;

[0196] 4) The position similarity score S pos , the mutation size similarity score S size , and the chromosome similarity score S Chr are weighted and summed to obtain the comprehensive similarity score S; expressed as:

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

[0198] Where w pos , w size are the weight coefficients of the position similarity score and the mutation size similarity score, w pos + w size = 1, for example 0.4, 0.6, the weight can be adjusted according to the specific analysis scene and data characteristics, so that the formula better adapts to the actual demand, and accurately measures the similarity between mutation characteristics.

[0199] II. For translocation mutation signals:

[0200]

[0201] Where S Chr1 represents the starting chromosome similarity score, represents the starting chromosome name of the first signal, represents the starting chromosome name of the second signal.

[0202]

[0203] wherein S Chr2 denotes the end chromosome similarity score, denotes the end chromosome name of the first signal, denotes the end chromosome name of the second signal;

[0204] the start chromosome similarity score S Chr1 , the end chromosome similarity score S Chr2 are weighted and summed to obtain a comprehensive similarity score S; denoted as:

[0205]

[0206] wherein, denotes the start position of the second signal, denotes the start position of the first signal, denotes the end position of the second signal denotes the end position of the first signal.

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

[0208] The other steps and parameters are the same as one of the first to sixth embodiments.

[0209] Embodiment eight: The embodiment is different from one of the first to seventh embodiments in that the step one five is performed to filter all the clustering cluster results outputted to obtain the final clustering results, i.e., the signals contained in each clustering cluster are integrated into one structural variation signal for output; the specific process is as follows:

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

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

[0213] 31) For the clustering results of the insertion variation type:

[0214] For each clustering cluster of the clustering results of the insertion variation type, if the number of clusters in the clustering cluster is less than or equal to max_cluster, calculate the average variation signal start position and average length of all variation data in the clustering cluster; the specific process is as follows:

[0215] The variation start positions of each variation signal in the first δ cluster are added one by one, and then divided by the number of variations to obtain the average variation start position; δ = 1, 2, …, N, N represents the total number of cluster clusters in the clustering result of the insertion variation type;

[0216] The lengths of all variations in the first δ cluster are added and summed, and then divided by the number of variations to calculate the average length;

[0217] The other information in the first δ cluster is kept as the same content as the first variation signal in the cluster cluster;

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

[0219] The average variation start position, the average length, and the other information are output as the first δ cluster result;

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

[0221] 32) For the clustering result of the deletion variation type:

[0222] For each cluster cluster of the clustering result of the deletion variation type, if the number of clusters in the cluster cluster is less than or equal to max_cluster, the average variation signal start position, the average variation signal end position, and the average length of all variation data in the cluster cluster are calculated; The specific process is:

[0223] The variation start positions of each variation signal in the first δ cluster are added one by one, and then divided by the number of variations to obtain the average variation start position; δ = 1, 2, …, N, N represents the total number of cluster clusters in the clustering result of the insertion variation type;

[0224] The variation end positions of each variation signal in the first δ cluster are added one by one, and then divided by the number of variations to obtain the average variation end position;

[0225] The lengths of all variations in the first δ cluster are added and summed, and then divided by the number of variations to calculate the average length;

[0226] The other information in the first δ cluster is kept as the same content as the first variation signal in the cluster cluster;

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

[0228] The average variation start position, the average variation end position, the average length, and the other information are output as the first δ cluster result;

[0229] until the output results of all cluster clusters of the clustering results of the deletion variation type are obtained, that is, the final deletion variation signal;

[0230] 33) Clustering results for inversion variation type:

[0231] For each cluster cluster of the clustering results of the inversion variation type, if the number of clusters in the cluster cluster is less than or equal to max cluster, the average variation signal start position and the average variation signal end position of all variation data in the cluster cluster are calculated; the specific process is:

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

[0233] The variation end positions of each variation signal in the δ"th cluster cluster are added one by one, and then divided by the number of variations to obtain the average variation end position;

[0234] The other information in the δ"th cluster cluster is kept as the same as the first variation signal in the cluster cluster;

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

[0236] The average variation start position, the average variation end position and the other information are taken as the output results of the δ"th cluster cluster;

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

[0238] 34) Clustering results for translocation variation type:

[0239] For each cluster cluster of the clustering results of the translocation variation type, if the number of clusters in the cluster cluster is less than or equal to max cluster, the average variation signal start position and the average variation signal end position of all variation data in the cluster cluster are calculated; the specific process is:

[0240] The variation start positions of each variation signal in the δ'"th cluster cluster are added one by one, and then divided by the number of variations to obtain the average variation start position; δ'" = 1, 2, …, N'", N'" represents the total number of cluster clusters in the clustering results of the translocation variation type;

[0241] The variation end positions of each variation signal in the δ'"th cluster cluster are added one by one, and then divided by the number of variations to obtain the average variation end position;

[0242] The other information in the δth cluster is kept as the same as the first variant signal in the cluster;

[0243] The other information is the start chromosome name, the end chromosome name, and the name from contig;

[0244] The average variant start position, the average variant end position, and the other information are output as the δth cluster result;

[0245] Until the output results of all the clusters of the translocation variant type clustering results are obtained, i.e., the final translocation variant signal;

[0246] 35), for the clustering results of the duplication variant type:

[0247] For each cluster of the clustering results of the duplication variant type, if the number of clusters in the cluster is less than or equal to max_cluster, the average variant signal start position and the average variant signal end position of all the variant data in the cluster are calculated; the specific process is as follows:

[0248] The variant start positions of the variant signals in the δth cluster are added one by one, and then divided by the number of variants to obtain the average variant start position; N represents the total number of clusters in the clustering results of the duplication variant type;

[0249] The variant end positions of the variant signals in the δth cluster are added one by one, and then divided by the number of variants to obtain the average variant end position;

[0250] The other information in the δth cluster is kept as the same as the first variant signal in the cluster; The other information is the chromosome name and the name from contig;

[0251] The average variant start position, the average variant end position, and the other information are output as the δth cluster result;

[0252]

[0253] Until the output results of all the clusters of the duplication variant type clustering results are obtained, i.e., the final duplication variant signal.

[0254] Because some of the signals are numerical types, such as the chromosome start position and the length, which can be averaged in the cluster, but part of the information, such as the chromosome name, is a string type, which cannot be averaged, so only the content kept as the same as the first variant signal in the cluster is needed, because these non-numerical types are consistent in the same cluster; ​​​​

[0255] The filtering mainly includes the following steps:

[0256] 1) In the post-processing stage of the clustering analysis, a strict number filtering operation is performed for each generated cluster. The default maximum cluster number is 5, which is used as a key criterion. After the clustering operation is completed and a series of cluster groups are obtained, the number of cluster groups is counted and judged.

[0257] 2) If the number of cluster groups obtained by counting exceeds the maximum cluster number, such cluster groups will be discarded and will not participate in the subsequent result generation process. The purpose of this is to avoid excessive and miscellaneous cluster groups introducing unnecessary noise data and interfering with the final analysis results, and to ensure the simplicity and reliability of the output results.

[0258] 3) Conversely, if the number of cluster groups is less than the default cluster number, a fine result extraction process is started. At this time, the average breakpoint position and the average length of all variant data in the cluster group are calculated respectively. Specifically, the breakpoint (breakpoint is the start position of the chromosome or the end position of the chromosome) positions of each variant in the cluster group are added one by one, and then divided by the number of variants, so as to accurately obtain the average breakpoint position; similarly, the lengths of all variants are added and divided by the number of variants to calculate the average length. Finally, the average breakpoint position and the average length calculated by the two calculations and the other information that remains the same as the first signal are used as the final output result corresponding to the cluster group.

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

[0260] In the clustering, the similarity calculation is based on a predefined feature similarity formula, which calculates the similarity score according to the properties of the features (such as position, size, chromosome); the threshold is a user-adjustable parameter used to define the similarity criterion between features. In actual operation, when the similarity score calculated according to the similarity formula between two features is higher than the threshold, it is determined that the two features have similarity and are classified into the same cluster group; otherwise, if the score is lower than the threshold, it is determined that the two features are quite different and are classified into different cluster groups. This method can automatically adapt to different sizes of feature sets, and has low computational complexity.

[0261] The other steps and parameters are the same as one of the first to seventh embodiments.

[0262] Specific implementation nine: different from one of the specific implementations one to eight, the step two is that the mixed sorting of the processed variation signal obtained in step one is performed, the spatial distribution clustering of the sorted variation signal is performed, it is judged whether the spatial distribution clustered variation signal is a potential complex variation, if yes, it is judged as a complex variation signal, the complex variation signal is obtained, and step three is performed; if not, it is judged as a simple variation signal, step three is performed, and the simple variation signal is one of insertion, deletion, inversion, translocation and duplication; the specific process is as follows:

[0263] For the complex variation part, the basic variation signal obtained in the last stage is mixed and sorted according to the position of the variation signal on the reference genome, so that the subsequent clustering analysis can be effectively grouped based on the spatial position; and then the clustering is performed according to the distance, and the potential complex variation site meeting the corresponding mode is further fused to obtain the final complex variation form.

[0264] 1) The processed variation signal obtained in step one is mixed and sorted according to the starting position on the reference genome in ascending order by using merge sorting or quick sorting;

[0265] After the first step, there is a container for each type of variation. Now, the several types need to be mixed together and sorted again according to the position on the reference genome. The method can choose merge sorting or quick sorting, but the sorted result is certain, that is, all the variations are arranged in ascending order according to the starting position of the variation on the reference genome.

[0266] The interval distance parameter Dist is introduced, and the default setting is 1000bp. This parameter is a key scale index for subsequent judgment of potential complex variation region.

[0267] 2) Taking the starting position of the first variation signal after sorting as the center of the genomic position, extending Dist / 2 base pairs from the center, and extending Dist / 2 base pairs from the center, a scanning window with a length of Dist is constructed;

[0268] All the variation signals in the scanning window are obtained;

[0269] From the starting of the first variation signal after sorting, the genomic position where the variation signal is located is taken as the center to construct a genomic interval with a length of Dist as a scanning window. Specifically, the window extends Dist / 2 base pairs forward and an equal amount of base pairs backward. Then, all the variation signals falling within the scanning window are sequentially checked.

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

[0271] The mutation types include insertion, deletion, inversion, translocation or duplication.

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

[0273] For each scanning window, the number of different types of structural mutation signals is strictly counted, and the structural mutation types include insertion, deletion, inversion, translocation and duplication. When the number of different types of structural mutation signals in a scanning window exceeds two, the region is determined as a potential complex mutation region.

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

[0275] The "CSV" is a common comma-separated value format identifier, which is convenient for data storage, transmission and subsequent analysis and processing.

[0276] Once a region is determined as a potential complex mutation region, the complex mutation mode fusion process is started. The system comprehensively collects the detailed information of all mutation signals in the region, including the mutation type, the chromosome name and the start position, and integrates a set of structured complex mutation description information, with the format defined as (CSV, mutation type, mutation start position). The "CSV" is a common comma-separated value format identifier, which is convenient for data storage, transmission and subsequent analysis and processing.

[0277] 5) Moving the scanning window by one mutation signal in sequence, and repeating 1)-4) until all complex mutation signals in the mixed sequence are processed.

[0278] After completing the scanning, determination, fusion and output process of the window corresponding to the first mutation signal, the scanning window is moved by one mutation signal in sequence, and the whole process is repeated seamlessly until all mutation signals in the sequence are processed.

[0279] The other steps and parameters are the same as one of the first to eighth embodiments.

[0280] The tenth embodiment is different from the first to ninth embodiments in that the mutation signal obtained in step 2 is genotyped in step 3, and the specific process is as follows:

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

[0282] Step three,

[0283] All the variation signal information obtained after the above steps includes chromosome name, variation start position, variation end position, variation length, and sequence information, etc., which varies slightly depending on different variation signal types.

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

[0285] For each type of variation signal, the initial variation signal is added to the clustering cluster, and the next variation signal is sequentially traversed until the clustering of the corresponding type of variation signal is completed:

[0286] For insertion (INS) type variation signals, three conditions are checked: one is that the chromosome name of the current variation signal and the chromosome name of the last variation signal in the current clustering cluster are the same; two is that the Euclidean distance between the chromosome start position of the current variation signal and the chromosome start position of the last variation signal in the current clustering cluster is less than 200bp; three is that the difference between the variation length of the current variation signal and the variation length of the last variation signal in the current clustering cluster is less than 50bp; if all three conditions are met, the current variation signal is added to the current clustering cluster; otherwise, a new clustering cluster is opened, and the current variation signal is added to the new clustering cluster;

[0287] For deletion (DEL), inversion (INV), or duplication (DUP) type variation signals, three conditions are checked: one is that the chromosome name of the current variation signal and the chromosome name of the last variation signal in the current clustering cluster are the same; two is that the Euclidean distance between the variation start position of the current variation signal and the variation start position of the last variation signal in the current clustering cluster is less than 200bp; three is that the Euclidean distance between the variation end position of the current variation signal and the variation end position of the last variation signal in the current clustering cluster is less than 200bp; if all three conditions are met, the current variation signal is added to the current clustering cluster; otherwise, a new clustering cluster is opened, and the current variation signal is added to the new clustering cluster;

[0288] For translocation (TRA) type of variant signal, four conditions are checked, one is that the chromosome name at the start position of the current variant signal is the same as the chromosome name at the start position of the last variant signal in the current cluster; two is that the chromosome name at the end position of the current variant signal is the same as the chromosome name at the end position of the last variant signal in the current cluster; three is that the Euclidean distance between the start position of the current variant signal and the start position of the last variant signal in the current cluster is less than 200bp; four is that the Euclidean distance between the end position of the current variant signal and the end position of the last variant signal in the current cluster is less than 200bp; if the four conditions are all met, the current variant signal is added to the current cluster; otherwise, a new cluster is started, and the current variant signal is added to the new cluster;

[0289] For complex variant (CSV) signal, three conditions are checked, one is that the chromosome name at the start position of the current variant signal is the same as the chromosome name at the start position of the last variant in the current cluster; two is that the Euclidean distance between the start position of the current variant signal and the start position of the last variant in the current cluster is less than 1000bp; three is that the type of the current variant signal is the same as the type of the last variant signal in the current cluster; if the three conditions are all met, the current signal is added to the current cluster; otherwise, a new cluster is started, and the current signal is added to the new cluster;

[0290] Through the above steps, each cluster is a distance cluster indicating potential same haplotype variants;

[0291] Step three two,

[0292] If the signal in the distance cluster is an insertion variant signal, the sequence edit similarity between the insertion variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm;

[0293] If the signal in the distance cluster is a deletion variant signal, the sequence edit similarity between the deletion variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm;

[0294] If the signal in the distance cluster is a duplication variant signal, the sequence edit similarity between the duplication variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm;

[0295] If the signal in the distance cluster is an inversion variant signal, the sequence edit similarity between the inversion variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm;

[0296] If the signal in the distance cluster is a translocation variant or a complex variant signal, step three four is directly executed;

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

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

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

[0300] For the duplication variant signal, the sequence between the variant start position and the variant end position of the variant signal in the reference genome with the same chromosome name as the duplication variant signal is extracted as the duplication variant sequence;

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

[0302] Step three,

[0303] According to the edit distance of any two signals in the obtained cluster, the signals in the cluster are operated as follows:

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

[0305] If the sequence edit similarity of the insertion variant sequences of the 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 one cluster; until all variant signals in the insertion cluster are judged;

[0306] If the sequence edit similarity of the deletion variant sequences of the 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 one cluster; until all variant signals in the deletion cluster are judged;

[0307] If the sequence edit similarity of the duplication variant sequences of the 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 one cluster; until all variant signals in the duplication cluster are judged;

[0308] If the sequence edit similarity of the inversion variant sequences of the 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 one cluster; until all variant signals in the inversion cluster are judged;

[0309] Step three four is executed;

[0310] Step three four, genotyping the variant signals in the clustering cluster; the specific process is:

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

[0312] If the number of variant signals in the clustering cluster is 1, it is determined that the heterozygous sequence, and the genotype is determined to be 0 / 1 or 1 / 0 according to the haplotype to which it belongs;

[0313] Otherwise, it is considered to be noise, and the corresponding clustering cluster is discarded;

[0314] Output the variant signal containing the genotyping information;

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

[0316] There are three genotypes, 0 / 1, 1 / 0, and 1 / 1; 1 represents that it is different from the reference genome and there is a variation; 0 represents 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 detection has two identical variations on the two haplotypes at this position; if it is 1 / 0, it means that the genome used for detection has only one variation on the first haplotype at this position; if it is 0 / 1, it means that the genome used for detection has only one variation on the second haplotype at this position.

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

[0318] The other steps and parameters are the same as one of the first to ninth embodiments.

[0319] The present application can also have other various embodiments, and those skilled in the art can make various corresponding changes and modifications according to the present application without departing from the spirit and essence of the present application, but these corresponding changes and modifications should all belong to the protection scope of the claims attached to the present application.

Claims

1. A method for variant calling based on end-to-end assembled genomes, characterized in that: The method specifically comprises the following steps: Step one, variant signal processing, to obtain processed variant signals; the specific process is: Step one, obtaining BAM format alignment information; Step one, filtering the BAM format alignment information obtained in step one, to obtain filtered BAM format alignment information; Step one, extracting the variant signals in the filtered BAM format alignment information; the specific process is: According to the flag flag in the filtered BAM format alignment information, the alignment information is divided into intra signals and inter signals; Signal processing is performed on the intra signals to obtain two kinds of insertion and deletion variant signals; Signal processing is performed on the inter signals to obtain five kinds of insertion, deletion, inversion, translocation and duplication variant signals; Step one, clustering the extracted variant signals, and outputting clustering cluster results; the specific process is: Unsupervised clustering and filtering method based on feature similarity: based on distance distribution, the sites close to each other and the features similar to each other are clustered together to form an average variant to achieve filtering effect; The specific process of the filtering is: 1) Set the maximum cluster number max_cluster; 2) If the cluster number in the cluster is greater than max_cluster, discard the corresponding cluster; Step two, mixed sorting of the processed variant signals obtained in step one, spatial distribution clustering of the sorted variant signals, judgment of whether the variant signals after spatial distribution clustering are potential complex variants, if yes, it is judged as a complex variant signal, and the complex variant signal is obtained, and step three is executed; if not, it is judged as a simple variant signal, and step three is executed; the simple variant signal is one of insertion, deletion, inversion, translocation and duplication; the specific process is: 1) Merge sorting or quick sorting algorithm is used to sort the processed variant signals obtained in step one in ascending order according to the starting position on the reference genome; 2) Taking the starting position of the first variant signal in the sorted variant signals as the center, extending Dist / 2 base pairs from the center, and extending Dist / 2 base pairs from the center, to form a scanning window with a length of Dist; All the variant signals in the scanning window are obtained; 3) judging the variant types of all the variant signals in the scanning window; The variant types include insertion, deletion, inversion, translocation or duplication; When the number of variant signal types in a scanning window exceeds two, it is determined that the scanning window is a potential complex variant region; otherwise, the scanning window is not a potential complex variant region; 4) Extracting detailed information of each complex variant signal in the potential complex variant region, the detailed information including the variant type, the variant chromosome name and the variant starting position of each variant signal, and the format of the complex variant signal is (CSV, chromosome name, variant type, variant starting position); Wherein, CSV is a format identifier; ​ 5), move the scanning window by one variant signal in sequence, repeat 1)-4) until all complex variant signals after mixed sequencing are completed traversal processing; Step three, genotype the variant signals obtained in step two; the specific process is: Step three, the different types of variant signals obtained in step two are respectively distance clustering; the specific process is: For the insertion type of variant signal, check three conditions: one is that the chromosome name of the current variant signal and the chromosome name of the variant in the current clustering cluster are the same; two is that the Euclidean distance between the starting position of the chromosome in the current variant signal and the starting position of the chromosome in the last variant signal in the current clustering cluster is less than 200bp; three is that the difference between the variant length of the current variant signal and the variant length of the last variant signal in the current clustering cluster is less than 50bp; if the three conditions are met, the current variant signal is added to the current clustering cluster; otherwise, a new clustering cluster is started, and the current variant signal is added to the new clustering cluster; For the deletion, inversion or repeat type of variant signal, check three conditions, one is that the chromosome name of the current variant signal and the chromosome name of the variant in the current clustering cluster are the same; two is that the Euclidean distance between the starting position of the current variant signal and the starting position of the last variant signal in the current clustering cluster is less than 200bp; three is that the Euclidean distance between the end position of the current variant signal and the end position of the last variant signal in the current clustering cluster is less than 200bp; if the three conditions are met, the current variant signal is added to the current clustering cluster; otherwise, a new clustering cluster is started, and the current variant signal is added to the new clustering cluster; For the translocation type of variant signal, check four conditions, one is that the chromosome name of the starting position of the current variant signal and the chromosome name of the starting position of the last variant signal in the current clustering cluster are the same; two is that the chromosome name of the ending position of the current variant signal and the chromosome name of the ending position of the last variant signal in the current clustering cluster are the same; three is that the Euclidean distance between the starting position of the current variant signal and the starting position of the last variant signal in the current clustering cluster is less than 200bp; four is that the Euclidean distance between the ending position of the current variant signal and the ending position of the last variant signal in the current clustering cluster is less than 200bp; if the four conditions are met, the current variant signal is added to the current clustering cluster; otherwise, a new clustering cluster is started, and the current variant signal is added to the new clustering cluster; For the complex variant signal, check three conditions, one is that the chromosome name of the starting position of the current variant signal and the chromosome name of the last variant starting position in the current clustering cluster are the same; two is that the Euclidean distance between the starting position of the current variant signal and the starting position of the last variant in the current clustering cluster is less than 1000bp; three is that the variant type of the current signal and the variant type of the last variant signal in the current clustering cluster are the same; if the three conditions are met, the current signal is added to the current clustering cluster; otherwise, a new clustering cluster is started, and the current signal is added to the new clustering cluster; Step three two, If the signals in the distance clustering cluster are insertion variant signals, the sequence edit similarity between the insertion variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm; If the signals in the distance clustering cluster are deletion variant signals, the sequence edit similarity between the deletion variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm; If the signals in the distance clustering cluster are repeat variant signals, the sequence edit similarity between the repeat variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm; If the signals in the distance clustering cluster are inversion variant signals, the sequence edit similarity between the inversion variant sequences of any two variant signals is calculated by using the minimum edit distance algorithm; If the signals in the distance clustering cluster are translocation variants or complex variants, steps three and four are directly executed; Step three, A similarity threshold is set; If the sequence edit similarity of the insertion variant sequences of the two variant signals is less than the similarity threshold, the two variant signals are classified into one clustering cluster; otherwise, the two variant signals are not one clustering cluster; until all variant signals in the insertion cluster are judged; If the sequence edit similarity of the deletion variant sequences of the two variant signals is less than the similarity threshold, the two variant signals are classified into one clustering cluster; otherwise, the two variant signals are not one clustering cluster; until all variant signals in the deletion cluster are judged; If the sequence edit similarity of the repeat variant sequences of the two variant signals is less than the similarity threshold, the two variant signals are classified into one clustering cluster; otherwise, the two variant signals are not one clustering cluster; until all variant signals in the repeat cluster are judged; If the sequence edit similarity of the inversion variant sequences of the two variant signals is less than the similarity threshold, the two variant signals are classified into one clustering cluster; otherwise, the two variant signals are not one clustering cluster; until all variant signals in the inversion cluster are judged; Step three and four are executed; Step three and four, genotype the variant signals in the clustering cluster; the specific process is: If the number of variant signals in the clustering cluster is 2 and the variant signals come from different haplotypes, it is determined that the same haplotype sequence is 1 / 1; If the number of variant signals in the clustering cluster is 1, it is determined that the heterozygous haplotype sequence is 0 / 1 or 1 / 0 according to the haplotype to which it belongs; Otherwise, it is considered to be noise, and the corresponding clustering cluster is discarded; Output the variant signals containing genotyping information; Write each variant signal record in VCF format into a file.

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

3. The end-to-end assembly genome-based variant detection method of claim 2, wherein: The step one two in the step one one obtains the BAM format comparison information is filtered, get filtered BAM format comparison information; The specific process is: According to the flag flag in the BAM format comparison information, the comparison information record is classified, and the type of the comparison information record is determined; The type of comparison information record can be divided into main comparison, secondary comparison or supplementary comparison; Through the SA label of main comparison, all fragment comparison information records constituting the same contig are obtained according to the same contig name, and formula 1 is given: contig = {seg1, seg2,..., seg j ,…,seg n}, 1 < i < n (1) The filtering condition is: wherein contig is alignment information of a contig, seg i is alignment information record of the i-th segment in contig sorted by relative position of segment start, ReadEnd segi represents the end position of seg i in contig; represents the end position of seg i-1 in contig; 1 < i ≤ n, n represents the total number of segments in contig; Selecting the alignment information as filtered BAM format alignment information.

4. The end-to-end assembly genome-based variant detection method of claim 3, wherein: The step one two in the step one one obtains the BAM format comparison information is filtered, get filtered BAM format comparison information; The specific process is: If the number in front of the I character in the CIGAR character bit in the filtered seg i 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. If the number in front of the D character in the CIGAR character bit in the filtered seg i is greater than 30bp, it is considered that there is a deletion signal, and the chromosome name where the deletion signal is located, the start position of the deletion signal on the reference genome, the end 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 end-to-end assembly genome-based variant detection method of claim 4, wherein: The step one two in the step one one obtains the BAM format comparison information is filtered, get filtered BAM format comparison information; The specific process is: Traverse the adjacent two segs constituting the contig, respectively recorded as seg1 and seg2, wherein the fragment with smaller starting coordinate is designated as the first aligned fragment; wherein represents the start position of seg1 on the stitching segment, represents the start position of seg2 on the stitching segment; denotes the end position of segl on the stitching segment, denotes the end position of seg2 on the stitching segment; denotes the start position of segl on the reference genome, denotes the start position of seg2 on the reference genome; denotes the end position of segl on the reference genome, denotes the end position of seg2 on the reference genome; read dis denotes the distance between the start position of seg2 on the stitched segment and the end position of seg1 on the stitched segment, ref dis denotes the distance between the start position of seg2 on the reference genome and the end position of seg1 on the reference genome, OLP ref denotes the overlap size on the reference genome between the end position of seg1 on the reference genome and the start position of seg2 on the reference genome; 1), for insertion variation signal, if the insertion variation signal satisfies the following precondition: wherein, represents the chromosome name in segl, represents the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in segl is consistent with the chromosome name recorded in seg2; And satisfy: Then extract the insertion signal: wherein Chrom represents the name of the chromosome where the extracted signal is located, Chrom INS represents the name of the chromosome where the extracted insertion signal is located, Start represents the start position identifier of the extracted signal, RefEnd seg1 represents the end position of the reference genome in seg1, RefEnd seg2 represents the end position of the reference genome in seg2, SVlen represents the length of the variation, Seq represents the sequence information identifier of the insertion signal, Seq INS represents the sequence corresponding to the insertion signal; If formula (4), (5) is not satisfied, the insertion signal is not extracted; 2), for deletion variation signal, if the deletion variation signal satisfies the following precondition: wherein, represents the chromosome name in segl, represents the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in segl is consistent with the chromosome name recorded in seg2; And satisfy: Then extract the deletion signal: wherein Chrom DEL represents the name of the chromosome where the extracted deletion signal is located, End represents the end position identifier of the extracted signal, and other identifiers are defined as above. If formula (4), (7) is not satisfied, the deletion signal is not extracted; 3), for repeat variation signal, if the repeat variation signal satisfies the following precondition: wherein, represents the chromosome name in segl, represents the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in segl is consistent with the chromosome name recorded in seg2; And satisfy: Then extract the repeat signal: wherein, ReadEnd seg1 denotes the end position of the splice junction in seg1, ReadStart seg2 denotes the start position of the splice junction in seg2, RefStart seg2 denotes the start position of the reference genome in seg2, Chrom DUP denotes the chromosome name where the extracted DUP signal is located, Th olp denotes the distance threshold for extracting the repetitive DUP signal; If formula (4), (9) is not satisfied, the repeat signal is not extracted; 4), for inversion variation signal; 41), if the inversion variation signal satisfies the following precondition: wherein, represents the chromosome name in segl, represents the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in segl is consistent with the chromosome name recorded in seg2; And satisfy: The alignment direction in seg1 and seg2 record is different, there is one reverse alignment and one forward alignment between seg1 and seg2, the forward aligned seg is recorded as fwd, and the reverse aligned seg is recorded as rev; | ReadEnd rev - ReadStart fwd | ≤ 100 bp (11) Then extract the inversion signal: If formula (4), (11) is not satisfied, the inversion signal is not extracted; 42), if the inversion variation signal satisfies the following precondition: wherein, represents the chromosome name in segl, represents the chromosome name in seg2; formula (4) indicates that the chromosome name recorded in segl is consistent with the chromosome name recorded in seg2; And satisfy: The alignment direction in seg1 and seg2 record is different, there is one reverse alignment and one forward alignment between seg1 and seg2, the forward aligned seg is recorded as fwd, and the reverse aligned seg is recorded as rev; | ReadEnd fwd - ReadStart rev | ≤ 100 bp (13) Extract the INV signal: If formula (4), (13) is not satisfied, the inversion signal is not extracted; wherein, ReadEnd rev is the start position of the reference genome in the forward alignment information, RefStart fwd is the start position of the reference genome in the reverse alignment information, RefStart rev is the start position of the reference genome in the forward alignment information, ReadEnd fwd is the start position of the reference genome in the reverse alignment information, ReadStart fwd is the start position of the reference genome in the forward alignment information, ReadStart rev is the start position of the reference genome in the reverse alignment information, RefEnd fwd is the start position of the reference genome in the forward alignment information, RefEnd rev is the start position of the reference genome in the reverse alignment information; Chrom INV is an identifier indicating the name of the chromosome where the INV variant signal occurs; Start INV is an identifier indicating the start position of the INV variant signal, End INV is an identifier indicating the end position of the INV variant signal; 5), for translocation variation signal; 51), if the translocation variation signal satisfies the following precondition: Wherein, 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 | ≤ 1000 bp (16) Then extract the translocation variation signal: If formula (15), (16) is not satisfied, the translocation signal is not extracted; 52), if the translocation variation signal satisfies the following precondition: Wherein, 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 | ≤ 1000 bp (18) Then extract the translocation signal: If formula (15), (18) is not satisfied, the translocation signal is not extracted; Wherein, Chrom1 represents the identification of the starting chromosome name of the TRA variant signal occurrence, Pos1 represents the identification of the starting position of the TRA variant signal occurrence, Chrom2 represents the identification of the ending chromosome name of the TRA variant signal occurrence, and Pos2 represents the identification of the ending position of the TRA variant signal occurrence.

6. The end-to-end assembly genome-based variant detection method of claim 5, wherein: The step one four clusters the extracted variant signals, and outputs a clustering cluster result; the specific process is as follows: (1) For the extracted deletion signal, insertion signal, duplication signal, inversion signal and translocation signal; Initialize the deletion clustering container; initialize the insertion clustering container; initialize the duplication clustering container; initialize the inversion clustering container; and initialize the translocation clustering container; (2) Sort the variant signals in each variant type set in ascending order according to the starting position in the reference genome; (3) For the signals in the deletion variant type signal set, the following operations are sequentially performed: a. Initialize an empty cluster, and add the first signal in the deletion variant type signal set as the starting signal to the cluster, and perform 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, add the current cluster to the deletion clustering container, and add the current signal to a new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal, and repeat b-c until all signals in the deletion variant type signal set are judged; (4) For the signals in the insertion variant type signal set, the following operations are sequentially performed: a. Initialize an empty cluster, and add the first signal in the insertion variant type signal set as the starting signal to the cluster, and perform 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, add the current cluster to the insertion clustering container, and add the current signal to a new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal, and repeat b-c until all signals in the insertion variant type signal set are judged; (5) For the signals in the duplication variant type signal set, the following operations are sequentially performed: a. Initialize an empty cluster, and add the first signal in the duplication variant type signal set as the starting signal to the cluster, and perform 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, add the current cluster to the duplication clustering container, and add the current signal to a new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal, and repeat b-c until all signals in the duplication variant type signal set are judged; (6) For the signals in the inversion variant type signal set, the following operations are sequentially performed: a. Initialize an empty cluster, and add the first signal in the inversion variant type signal set as the starting signal to the cluster, and perform 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, add the current cluster to the inversion clustering container, and add the current signal to a new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal, and repeat b-c until all signals in the inversion variant type signal set are judged; (7) For the signals in the translocation variant type signal set, the following operations are sequentially performed: a. Initialize an empty cluster, and add the first signal in the translocation variant type signal set as the starting signal to the cluster, and perform 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, add the current cluster to the translocation clustering container, and add the current signal to a new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal, and repeat b-c until all signals in the translocation variant type signal set are judged. 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, add the current cluster to the inversion cluster container and add the current signal to the new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal and repeat b-c until all signals in the set of variant signals of the inversion mutation type are judged; (7) For the signals in the set of variant signals of the translocation mutation type, sequentially perform the following operations: a. Initialize an empty cluster and add the first signal in the set of variant signals of the translocation mutation type as the starting signal, and perform 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, add the current cluster to the translocation cluster container and add the current signal to the new cluster; If the similarity is greater than or equal to the predetermined threshold, add the current signal to the current cluster; d. Move to the next signal and repeat b-c until all signals in the set of variant signals of the translocation mutation type are judged; (8) Obtain all cluster results of all mutation type signals.

7. The end-to-end assembly genome-based variant detection method of claim 6, wherein: The similarity between each signal and the last signal of the current cluster is calculated in b; the specific process is: I. For insertion, deletion, inversion, and duplication variant signals: 1) Calculate position similarity S pos score; the specific process is: Position similarity S pos The score is calculated as follows: S pos = |start1 - start2| (20) Where start1 and Start2 are the starting positions of the two variant signals; 2) Calculate the variant signal size similarity score S size ; the specific process is: The variant length sizes of the two variant signals are SVlenl and SVlen2, respectively, and the variant length size similarity score S of the variant signals is calculated by the following formula: size is calculated by the following formula: 3) Calculate the chromosomal similarity score S Chr ; the specific process is: Let the two chromosomes where the variant signals are located be Chroml and Chrom2, respectively, and the chromosome similarity score S Chr is calculated by the following discrete function: If the chromosome names of the two variant signals are the same, the chromosome similarity score is 1; If the chromosome names of the two variant signals are different, the chromosome similarity score is 0; 4) the position similarity S pos score, the variant signal size similarity score S size score, the chromosome similarity score S Chr are weighted and summed to produce a composite similarity score S; denoted as: S = S Chr x (w size x S size + w pos x S pos ) (23) wherein w pos , w size are weight coefficients of the position similarity score, the variant signal size similarity score, respectively, w pos +w size = 1. II. For translocation variant signals: where S Chr1 represents the starting chromosome similarity score, represents the starting chromosome name of the first signal, represents the starting chromosome name of the second signal; where S Chr2 denotes the end chromosome similarity score, denotes the end chromosome name of the first signal, denotes the end chromosome name of the second signal; The starting chromosomal similarity score S Chr1 The ending chromosomal similarity score S Chr2 The weighted sum is taken to arrive at a composite similarity score S; denoted as: wherein denotes the start position of the second signal, denotes the start position of the first signal, denotes the end position of the second signal, denotes the end position of the first signal; w pos1 is a weight factor for the start position, w pos2 is a weight factor for the end position.

8. The end-to-end assembly genome-based variant detection method of claim 7, wherein: The step one five filters all the output cluster results to obtain the final clustering result, that is, the signals in each cluster are integrated into a structural variant signal for output; the specific process is: 1) Set the maximum cluster number max_cluster; 2) If the number of clusters in the cluster is greater than max_cluster, discard the corresponding cluster; 3)、 31) For the clustering results of the insertion mutation type: For each cluster in the clustering results of the insertion mutation type, if the number of clusters in the cluster is less than or equal to max_cluster, calculate the average variant signal starting position and average length of all variant data in the cluster; the specific process is: Add the variant starting positions of the variant signals in the δth cluster one by one, and then divide by the number of variants to obtain the average variant starting position; δ = 1, 2, …, N, N represents the total number of clusters in the clustering results of the insertion mutation type; Sum the lengths of all variants in the δth cluster and divide by the number of variants to calculate the average length; Keep the other information in the δth cluster the same as the first variant signal in the cluster; The other information is the chromosome name, the name from the contig, and the inserted sequence; Output the average variant starting position, average length, and other information as the δth cluster output result; Until the output results of all clusters of the clustering results of the insertion mutation type are obtained, that is, the final insertion variant signal; 32) For the clustering results of the deletion mutation type: For each cluster of the clustering result of the deletion mutation type, if the number of clusters in the cluster is less than or equal to max cluster, the average mutation signal start position, the average mutation signal end position and the average length of all the mutation data in the cluster are calculated; the specific process is as follows: The mutation start positions of all the mutation signals in the δ'th cluster are added 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 the deletion mutation type; The mutation end positions of all the mutation signals in the δ'th cluster are added one by one, and then divided by the number of mutations to obtain the average mutation end position; The lengths of all the mutations in the δ'th cluster are added and divided by the number of mutations to obtain the average length; The other information in the δ'th cluster is kept as the same as that of the first mutation signal in the cluster; The other information is the name of the chromosome and the name from the contig; The average mutation start position, the average mutation end position, the average length and the other information are output as the δ'th cluster; Until the output results of all the clusters of the clustering result of the deletion mutation type are obtained, that is, the final deletion mutation signal; 33) For the clustering result of the inversion mutation type: For each cluster of the clustering result of the inversion mutation type, if the number of clusters in the cluster is less than or equal to max cluster, the average mutation signal start position and the average mutation signal end position of all the mutation data in the cluster are calculated; the specific process is as follows: The mutation start positions of all the mutation signals in the δ''th cluster are added 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 the inversion mutation type; The mutation end positions of all the mutation signals in the δ''th cluster are added one by one, and then divided by the number of mutations to obtain the average mutation end position; The other information in the δ''th cluster is kept as the same as that of the first mutation signal in the cluster; The other information is the name of the chromosome and the name from the contig; The average mutation start position, the average mutation end position and the other information are output as the δ''th cluster; Until the output results of all the clusters of the clustering result of the inversion mutation type are obtained, that is, the final inversion mutation signal; 34) For the clustering result of the translocation mutation type: For each cluster of the clustering result of the translocation mutation type, if the number of clusters in the cluster is less than or equal to max cluster, the average mutation signal start position and the average mutation signal end position of all the mutation data in the cluster are calculated; the specific process is as follows: The mutation start positions of all the mutation signals in the δ'''th cluster are added 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 the translocation mutation type; The mutation end positions of all the mutation signals in the δ'''th cluster are added one by one, and then divided by the number of mutations to obtain the average mutation end position; Keeping other information in the δ'''th cluster as the same as the first variant signal in the cluster; The other information is the start chromosome name, the end chromosome name, and the name from contig; Output the average variant start position, the average variant end position, and the other information as the output result of the δ'''th cluster; Until the output result of all clusters of the clustering result of the translocation variant type is obtained, that is, the final translocation variant signal; 35), for the clustering result of the duplication variant type: For each cluster of the clustering result of the duplication variant type, if the number of clusters in the cluster is less than or equal to max_cluster, the average variant signal start position and the average variant signal end position of all variant data in the cluster are calculated; the specific process is: The variant start positions of each variant signal in the first cluster are accumulated one by one, and then divided by the number of variants to obtain the average variant start position; The variant start positions of each variant signal in the first cluster are accumulated one by one, and then divided by the number of variants to obtain the average variant start position; N"'represents the total number of cluster clusters in the clustering result of the repeated variant type; The termination positions of the variants in the first cluster are added one by one, and then divided by the number of variants to obtain the average termination position of the variants. of the variants in the first cluster are added one by one, and then divided by the number of variants to obtain the average termination position of the variants. The other information in the first cluster is kept as the first variant signal in the cluster. ​ The other information is the chromosome name and the name from contig; The average variant start position, average variant end position, and other information are output as the first cluster output result, Until the output result of all clusters of the clustering result of the duplication variant type is obtained, that is, the final duplication variant signal.

Images